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Matched  Field  Processing  -  What  Next? 


John  A.  Tague 

Undersea  Surveillance  Signal  Processing  Team 
Office  of  Naval  Research 
800  N.  Quincy  St. 

Arlington,  VA  22217-5660,  USA 

12  June  1996 


Abstract 

Matched  field  processing  is  an  intellectually  fascinating  blend  of  statistics  and  underwater  acoustics.  Further¬ 
more,  it  works!  Many  experiments  have  demonstrated  its  viability  -  in  both  shallow  water  and  deep  ocean 
environments.  Where  do  we  go  next?  What  challenging  problems  remain  to  be  solved,  and  what  options 
can  we  offer  the  U.S.  Navy  for  practical  arrays  capable  of  supporting  this  advanced  processing  technique? 
The  Office  of  Naval  Research  is  interested  in  developing  matched  field  processing  that  can  provide  a  new 
and  vital  fieet  capability. 
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limiting  factors:  noise,  reverberation 


one  way  r(t)  =  h(t)*s(t)  +  n(t) 


Signal  Detection 
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1 .  correlate  s,  r  {matched-filtet)  for  all  lags 

2.  obtain  the  maximum  of  the  envelope  of  the 
correlation 


Signal  Detection 
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Whalen,  Detection  of  Signals  in  Noise 
Johnson  and  Dudgeon,  Array  Signal  Processing 


Signal  Detection 
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(a)  model  based  matched  filter 

(b)  model  based  matched  filter  with  uncertain 
amplitude 

(c)  standard  matched  filter 
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3.  Is  the  pay  off  from  the  computation  of  the 
optimum  processor  significant? 


Simulations  Scenario 
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Frequency:  300  -  900  Hz 
Background:  white  Gaussian  noise 


Dispersion  Curves  of 
Waveguide 
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1540  1530  1520  1510  1500  1490 

group  velocity  (m/s) 


Dispersion 


time-series  for  1 2  km  range 

sine  pulse  (  500  -  900  Hz ) 

We  can  identify  pulses  carried 
through  distinct  modes. 


500  -  900  Hz,  12  km 
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Dispersion  -  Distortion 
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depth  (m)  depth  (m) 


time  (s) 

sr  =  2  km,  modes  have  not  separated 

yet- 
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crcBs-correiaiiOT  (standard  matched-filter) 
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frequency  ,  less  dispersion  low  frequency  ,  more  dispersion 


depth  (m)  depth  (m) 


Correlation  vs.  Frequency 

(standard  matched-filter) 


500  -  900  Hz 
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Ambiguity  Surfaces  of  Standard 

Matched-Fiiter 


ROC  Performance  Evaluation 
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d  =  Es  /  a2 

Es:  signal  energy,  a2=noise  variance 


■ 


CO 

> 

LU 

0) 

O 


o 

o 

oc 


0>C0N<D10'^C0CM'i-O 

ooddddddd 


00 

d 


CO 

d 


d 


cvi 

d 


o 


uoipa^ep  |o  X^inqBqojd 


(uj)  mdap 


21 


time  (sec)  probability  of  false  alarm 

standard  matched-filter 


Detection  Performance 
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probability  of  false  alarm  probability  of  false  alarm 


Detection  Performance 


probability  of  false  alarm  probability  of  false  alarm 
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Optimum  MBMF  Detection 
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uenetic  Aigoritnms  (p.  Gerstoft) 

Simulated  Annealing  (Collins,  Kuperman  -  N.  Frazer 
Michalopoulou,  Porter) 


Summary 
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requires  accurate  modeling  of  the  involved  uncertainty 
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ACTIVE  MATCHED-FIELD  TRACKING  (AMFT) 


Homer  Bucker,  NRaD,  San  Diego  CA 


A  primary  objective  of  a  shallow  water  ASW  system  is  to  develop  tracks  of  an  under¬ 
water  submarine.  In  matched-field  tracking  (J.  Acoust.  Soc.  Amer.  96, 3809-,  1994) 
outputs  of  a  set  of  passive  sensors  are  processed  to  find  the  best  match,  i.e.  highest 
correlation,  between  elements  of  the  measured  covariance  matrix  and  those  calculated  for 
possible  tracks.  The  processing  is  straight-forward  but  computer  intensive  because  all 
data  collected  over  a  time  window  of  several  minutes  and  a  frequency  band  of  several 
hundred  Hz  must  be  carefully  examined. 

In  the  active  case,  similar  methods  can  be  used.  However,  the  definition  of  a 
covariance  matrix  needs  to  be  extended  (J.  Acoust.  Soc.  Amer.  99, 1784-,  1996). 
Fortunately,  only  a  small  fraction  of  the  data  collected  need  be  compared  to  calculated 
values.  This  is  because  with  a  given  set  of  pulses  and  a  postulated  track,  only  a  few 
received  values  are  compared  with  those  calculated. 

The  matched-field  tracking  algorithm  was  tested  with  a  simulation  of  an  ASW  system 
operating  in  a  200  m  shallow  water  channel  with  a  sloping  bottom.  As  shown  in  the 
figure,  three  surface  ships  and  a  submarine  at  100  m  depth  were  pinged  upon  (160  pings) 
for  a  period  of  5  minutes.  The  source  was  a  vertical  string  of  10  transducers  emitting  500 
ms  CW  pings  at  a  frequency  of  250  Hz  and  was  located  at  the  origin  of  the  coordinate 
system.  When  transmitting,  each  transducer  radiated  1 0  w  of  power  so  that  the  average 
radiated  power  in  the  water  is  less  that  3  watts.  The  receiving  array  consisted  of  7  sensors 
collocated  with  the  transducer  string.  There  were  7  hydrophones,  with  5  sensors  located 
along  the  X  axis  and  2  sensors  along  the  Y  axis.  The  ambient  noise  level  weis  65  dB/pPa. 

The  results  of  the  simulation  are  shown  in  the  figure  where  the  best  track  for  an 
xmder-water  target  is  shown  as  the  curve  of  dots  that  closely  follows  the  submarine  track 
(the  solid  curve).  The  estimated  track  is  a  composite  of  16  associated  segments.  Each 
segment  is  defined  by  six  coordinates,  the  x  and  y  values  of  the  end  points,  a  constant 
depth,  and  a  measure  of  the  curvature.  This  track  has  a  correlation  value  of  0.178  and  the 
estimated  target  depth  varied  between  98  and  107  m.  The  simulated  depth  was  100  m. 

In  the  current  example,  10  multi-paths  were  used  to  calculate  the  simulated  signals 
received  at  the  array  and  2  multi-paths  were  used  to  calculate  the  sound  pressures  for  the 
covariance  elements  associated  with  possible  tracks  While  the  results  seem  promising,  at 
sea  experiments  are  needed  to  validate  the  AMFT  concept. 
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Matched-Field  Track-Before-Detect  (TBD)  Processing 

using  SWellEX  Data 

Paul  A.  Baxley 

Naval  Command  Control  and  Ocean  Surveillance  Center 
Research,  Development,  Test  and  Evaluation  Division 
San  Diego,  CA  91935  USA 

12  June  1996 


Abstract 

Matched-field  tracking  is  a  generalization  of  matched-field  processing  in  which  the  time  evolution  of 
matched-field  ambiguity  surfaces  are  used  in  a  source-track  search.  Track-before-detect  (TBD)  processing 
makes  use  of  this  technique  to  extract  source  track  information  so  that  the  source  signal  function  can  be 
reconstructed,  rendering  it  more  detectable.  The  detection  enhancement  produced  by  a  shift-then-average 
TBD  algorithm  is  investigated  for  multitone  source-tow  data  recorded  near  San  Diego  during  the  first 
and  third  shallow  water  evaluation  cell  experiment  (SWellEX-1  and  3).  The  source- tow  (approximately 
3  knots  in  SWellEX-1  and  5  knots  in  SWellEX-3)  traversed  a  range-independent  radial  from  a  FLIP- 
mounted  vertical  line  array  (VLA).  Previously  masked  sources  were  rendered  detectable  in  ambient  noise 
(SWellEX-3)  and  in  the  presence  of  a  single  loud  interferer  (SWellEX-1).  Frequency  averaging  appears 
to  improve  TBD  performance  primarily  for  cases  in  which  the  source  peak  is  masked  by  strong  sidelobes 
of  noise  sources,  as  occurs  when  an  interferer  passes  nearby. 


Introduction 

Matched-field  tracking  [1]  is  a  generalization  of  matched-field  processing  in  which  the  time  evolution  of 
matched-field  ambiguity  surfaces  are  used  in  a  source-track  search.  Matched-field  track-before-detect  (TBD) 
processing  [2], [3]  makes  use  of  this  technique  to  extract  source  track  information  so  that  the  source  signal 
function  can  be  reconstructed,  rendering  it  more  detectable.  For  a  vertical  line  array  in  a  range-independent 
environment  and  a  constant  range-rate  source,  this  may  be  accomplished  via  a  simple  shift-then-average 
approach,  which  3delds  the  track  parameters  corresponding  to  the  optimal  source  track.  The  objective  of 
this  work  is  to  experimentally  demonstrate  the  feasibility  of  this  approach  and  to  evaluate  the  performance 
gains  which  may  be  achieved.  This  was  accomplished  via  the  analysis  of  multitone  source-tow  data  recorded 
near  San  Diego  during  the  first  and  third  shallow  water  evaluation  cell  experiment  (SWellEX-1  and  3). 


I  Scenario  assumptions 

For  the  most  general  case  of  an  asymmetric  environment  or  array,  matched-field  trackers  must  search  for 
optimal  tracks  from  a  set  of  those  defined  by  curves  connecting  all  combinations  of  potential  start  and  end 
points.  If  an  assumption  of  constant  speed  is  assumed,  these  curves  collapse  to  straight  lines.  A  maximum 
correlation  between  measured  and  predicted  pressures  arriving  at  the  array  sensors  over  a  given  time  interval 
yields  the  optimal  track,  characterized  by  6  track  parameters;  namely,  the  three  cartesian  coordinates  of  the 
start  and  end  points  of  the  track.  A  full  search  for  this  case  is  obviously  a  formidable  problem  and  will  often 
be  computationally  restrictive. 
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A  significant  reduction  of  the  search  space  occurs  for  the  case  of  VLA  processing  in  a  symmetric  (range- 
independent)  environment.  For  such  a  scenario,  the  array  cannot  distinguish  the  azimuthal  direction  of 
incoming  energy.  Consequently  the  search  is  transformed  into  a  search  for  a  family  of  tracks,  which  may 
occur  at  any  azimuth.  The  search  space  for  this  Ccise  is  reduced  to  four  parameters:  the  range  and  depth 
at  either  the  start  or  end  of  the  track,  the  speed,  and  the  range  at  the  closest  point  of  approach  (CPA). 
If  the  search  is  further  restricted  to  tracks  along  a  radial  line  connecting  the  source  and  VLA,  which  is 
equivalent  to  a  constant  range-rate  assumption,  then  the  range  at  CPA  is  zero  and  the  track  parameters  are 
reduced  to  three:  the  range  and  depth  at  either  the  start  and  end  of  the  track,  and  the  range-rate.  This  last 
scenario  was  selected  for  this  analysis  for  several  reasons.  First,  the  computational  load  is  greatly  reduced 
allowing  for  more  efficient  evaluation  of  the  TBD  techniques.  Second,  good  quality  SWellEX  data  with  such 
scenarios  is  readily  available.  In  addition,  the  SWellEX  tracks  analysed  in  this  study  are  along  nearly  constant 
bathymetry  lines,  so  that  range-dependent  or  three-dimensional  effects  need  not  be  considered.  Third,  this 
scenario  allows  for  the  use  of  a  simple  shift-then-average  TBD  technique,  to  be  discussed  shortly.  And 
lastly,  w'hile  this  shift-then-average  technique  implies  on-radial  tracks,  it  can  still  be  applied  approximately 
to  off-radial  tracks  at  long  range,  since  the  range-rate  varies  slowly  in  this  case. 


II  Processing  Technique 


Given  the  above  track  scenario  assumption,  the  optimal  track  can  be  obtained  by  a  simple  shift- then-average 
process  applied  to  consecutive  matched-field  ambiguity-surface  time  samples.  Assume  that  over  a  given  time 
interval  M  samples  of  acoustic  pressure  are  measured  at  each  phone  of  the  VLA.  A  measure  of  the  match 
between  the  observed  pressure  at  each  phone  of  an  iV-phone  array  and  the  predicted  pressure  Pn(x)  at 
the  same  phones  for  an  assumed  source  location  vector  x  is  given  by  the  Bartlett  estimator,  defined  as 


Bbart{^)  = 


En=l  Em=l  Pm(x)iUnPn(x) 

Etl|Pn(x)pE[EtiKr]' 


(1) 


where  E[]  denotes  the  time  average  (expected  value),  and  Rmn  =  E\p^p^*],  the  time-averaged  cross- 
spectral  matrix  of  the  observed  pressures.  Application  of  (1)  to  each  measured  time  sample  yields  a  range- 
depth  ambiguity  surface  for  each  sample.  For  a  source  moving  at  a  constant  depth,  the  range  of  the  source 
peak  will  change,  according  to  its  motion  trajectory,  for  each  sample.  Averaging  the  surfaces  without  taking 
into  account  this  motion  will  therefore  cause  a  degradation  in  the  correlation  level  of  the  source  peak.  If 
however,  the  surfaces  are  shifted  in  the  range  direction  by  the  amount  the  source  has  moved  for  each  time 
interval,  the  source  peaks  in  each  surface  will  be  aligned  and  the  source  peak  will  experience  no  degradation 
upon  averaging.  By  searching  through  different  range  rates,  which  de&es  the  shift  between  surfaces,  the 
range-rate  at  which  a  maximum  average  correlation  occurs  can  be  determined.  In  this  way  the  gain  in 
detectability  achievable  by  averaging  for  a  stationary  source  (blogM)  can  be  regained  for  the  case  of  the 
moving  source. 


Ill  SWellEX  1  and  3  environment 

SWellEX  1  and  3  took  place  near  the  southern  California  coast  in  approximately  50-200  m  water  southwest 
of  Point  Loma,  the  entrance  to  San  Diego  Harbour.  The  data  analyzed  in  this  work  is  that  measured  on 
Marine  Physical  Laboratory  (MPL)  SRP  VLAs  during  mutlitone  source  tows  along  tracks  of  nearly  constant 
water  depth  (198  m)  and  flat  bath}unetry.  Figure  1  shows  the  experimental  area  along  with  the  VLA  location 
and  source  tracks  performed  on  a  single  day  (August  17,  1993)  during  SWellEX-1.  The  VLA  location  during 
SWellEX-3  was  nearly  identical  to  that  in  SWellEX-1.  Only  the  track  northward  from  the  VLA  in  Fig.l, 
and  a  similar  northward  track  during  SWellEX-3,  are  considered  in  this  work.  For  each  data  set,  80  minutes 
of  data  was  processed  using  8192-point  FFT’s  with  50%  overlap,  resulting  in  a  total  of  1775  FFT  samples. 
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Depth 

(m) 

Compres- 

sional 

Wave 

Speed 

(m/s) 

Shear 

Wave 

Speed 

(m/s) 

Density 

(g/cm®) 

Compres- 

sional 

Wave 

Attenuation 

(dB/(kmHz)) 

Shear 

Wave 

Attenuation 

(dB/(kmHz)) 

Bottom 

Type 

198.0 

1572.368 

0.0 

1.76 

0.20 

0.0 

sediment 

228.0 

1593.016 

0.0 

1.76 

0.20 

0.0 

228.0 

1881.000 

0.0 

2.06 

0.06 

0.0 

mudstone 

1028.0 

3245.800 

0.0 

2.06 

0.06 

0.0 

1028.0 

5200.000 

0.0 

2.66 

0.02 

0.0 

— 1 

basement 

Table  1:  Geoacoustic  Model  for  SWeUEXl,3. 


With  a  sampling  rate  of  1500  sample/sec,  the  FFT  binwidth  was  0.1831  Hz,  the  FFT  length  was  5,46  sec, 
and  the  time  between  consecutive  samples  was  2.73  sec. 

The  SRP  VLA  contained  48  elements  with  a  88.125-m  aperture  for  SWellEX-l  and  63  elements  with  a 
116.25-m  aperture  for  SWellEX-3.  The  element  spacing  was  1.875  m  in  both  cases.  The  SWellEX-1  analysis 
considered  the  full  48-element  VLA,  whereas  the  SWellEX-3  analysis  used  only  a  subset  of  8  elements,  spaced 
by  11.25  m,  with  a  total  aperture  of  90  m.  The  reason  for  using  the  sparser  array  in  the  latter  case  was 
the  discovery  in  previous  SWellEX-1  analyses  [4]  that  there  are  only  about  17  propagating  modes  present  in 
the  data  at  the  highest  frequency  (195  Hz),  thus  requiring  a  sparser  vertical  sampling  of  the  field.  While  a 
number  of  tonals  were  transmitted  in  each  experiment,  only  a  strong  95-Hz  tonal  in  SWellEX-1  and  a  weak 
125-Hz  tonal  in  SWellEX-3  are  considered  in  this  analysis. 

Superimposed  on  Fig.  1  is  the  sound  speed  profile  measured  near  the  VLA  site.  Measurements  at  other 
sites  and  times  varied  little  firom  this  profile.  The  geoacoustic  model  in  Table  1  is  thought  to  adequately 
represent  the  bottom  properties  for  this  region  and  was  used  in  the  generation  of  all  replicas.  This  model  is 
close  to  that  determined  by  inversion  studies  to  be  optimal  [4]. 


IV  Experimental  results 

The  95- Hz  tonal  in  the  SWellEX-1  data  was  selected  for  analysis  primarily  because  the  occurrence  of  periods 
of  very  high  signal-to-noise  ratio,  which  facilitate  an  assessment  of  the  gain  in  detectibility  obtainable  via 
the  TBD  technique.  However,  the  presence  of  a  strong  unknown  interferer  during  a  portion  of  the  track  also 
allowed  for  a  performance  analysis  of  the  practical  case  of  a  low  signal-to-interferer-ratio  scenario.  The  125- 
Hz  tonal  in  the  SWellEX-3  data  was  one  of  the  quietest  tonals  of  that  experiment,  providing  for  a  practical 
analysis  of  the  TBD  technique  at  low  signal-to-noise  ratios  (between  -2  and  -12  dB).  The  shift-then-average 
TBD  technique  was  applied  to  these  80-min  range-independent  data  sets  using  averaging  times  of  0.32  min 
(8  samples),  0.68  min  (16  samples),  1.41  min  (32  samples),  2.87  min  (64  samples),  5.78  min  (128  samples). 
Averaging  intervals  were  overlapped  by  75%  in  all  cnses.  The  SWellEX-1  results  will  be  presented  first, 
followed  by  the  SWellEX-3  results. 

IV.l  SWellEX-1  results 

The  TBD  technique  searches  through  potential  range-rates,  via  the  shift-then-average  process,  to  obtain 
the  range-rate  and  source  range  and  depth  at  the  start  or  end  of  the  optimal  track  corresponding  to  the 
maximum  correlation  in  the  set  of  averaged  surfaces.  The  results  of  the  range-rate  search  for  the  SWellEX-1 
data  are  summarized  in  Fig.  2,  which  plots  maximum  correlation  in  an  averaged  surface  (corresponding  to 
a  particular  range  rate)  versus  the  number  of  ambiguity  surface  samples  included  in  the  average  for  the  full 
duration  of  the  source-tow  track.  Figure  2a  corresponds  to  16  samples  (0.68  min),  Fig.  2b  to  32  samples 
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(1.41  min),  Fig.  2c  to  64  samples  (2.87  min),  and  Fig.  2d  to  128  samples  (5.78  min).  The  first  column  is 
for  the  95-Hz  tonal,  while  the  second  column  is  for  a  noise-only  bin  centered  at  96  Hz,  included  to  facilitate 
an  understanding  of  tracker  in  this  noise  environment.  Each  point  in  these  plots  has  a  maximum  correlation 
level,  range-rate,  and  source  range  and  depth  (at  the  end  of  potential  track)  associated  with  it.  It  is  first 
observed  for  the  95-Hz  tonal  that  with  a  low  number  of  samples  (8  to  32)  a  track  is  clearly  discernible  for 
times  less  than  25  min  and  greater  than  about  45  min.  This  gap  in  track  detection  is  caused  by  the  presence 
of  an  interferer  during  this  period.  The  noise-only  plots  (second  column)  clearly  indicate  the  presence  of  a 
high-speed  (>  10  m/s)  track  in  the  noise,  which  masks  the  source-tow  track  in  the  signal+noise  plots  (first 
column).  In  addition,  even  though  the  track  is  detectable  for  the  low  number  of  averages,  the  range-rate 
resolution  is  poor.  Increasing  the  averaging  time  (number  of  samples)  however,  simultaneously  increases  the 
range-rate  resolution  and  renders  the  track  more  detectable  during  the  time  of  the  interferer.  These  effects 
are  more  clearly  observed  in  the  maximum  correlation  versus  range-rate  curves  obtained  at  specific  times 
in  Fig.  3.  Figure  3a  plots  the  results  for  the  case  in  which  the  averaging  time  begins  at  30  min  from  the 
start  of  the  track,  during  the  period  in  which  the  interferer  is  present,  while  Fig.  3b  plots  the  results  for 
the  case  in  which  the  averaging  time  begins  at  52  min  from  the  start  of  track,  when  the  interferer  is  absent 
and  signal-to-noise  level  is  high.  When  the  interferer  is  present  and  the  averaging  time  is  small,  peaks  in 
the  maximum  correlation  versus  range-rate  curves  tend  toward  the  higher  range  rates,  under  the  influence  of 
the  interferer,  as  observed  in  the  noise-only  curves.  However,  when  the  number  of  samples  exceeds  32  (1.41 
min),  the  source  peaks  along  the  true  track  have  reinforced  each  other  to  the  point  of  being  detectable  above 
the  interferer’s  track,  resulting  in  a  peak  at  the  true  range  rate  of  1.6  m/s.  The  maximum  correlation  versus 
range  rate  curves  for  the  high  signal-to-noise  period  following  the  period  of  the  interferer  demonstrates  nicely 
the  improved  range-rate  resolution  with  increased  averaging  time.  The  noise-only  results  during  this  period 
indicate  no  preferential  range  rate,  except  at  the  higher  averaging  times  when  the  peaks  tend  toward  a  zero 
range  rate.  This  tendency  towards  a  zero  range  rate  in  the  noise  is  also  observed  in  Fig.  2,  for  periods  before 
and  after  the  interferer.  This  effect  is  probably  caused  by  the  fact  that  source  of  noise  during  these  periods 
is  long-range  shipping,  whose  range-rates  vary  slowly  over  the  averaging  times  considered.  However,  as  a 
ship  moves  into  short  range,  as  with  the  interferer,  peaks  move  away  from  the  zero  range-rate  position. 

In  practice,  the  tracker  could  be  automated  to  report  the  maximum  correlation  from  among  the  set  of 
averaged  surfaces  (corresponding  to  candidate  range-rates)  and  the  range-rate,  source  range,  and  source 
depth  associated  with  that  maximum.  Figure  4  presents  such  automated  results,  as  a  function  of  number  of 
samples,  over  the  source-tow  period.  Each  dot  in  this  figure  represents  the  result  for  an  averaging  interval, 
while  the  solid  curves  in  the  range  and  depth  versus  time  plots  represent  the  known  values.  Note  that  for 
the  small  averaging  times,  while  there  is  considerable  scatter,  there  is  a  definite  congregating  of  data  near 
the  true  values,  except  when  the  interferer  is  present  (between  25  and  45  min).  Also,  note  that  the  speed 
(or  range  rate)  data  is,  on  average,  offset  from  the  zero  range-rate  position  when  the  interferer  is  absent,  but 
is  scattered  over  a  wide  range  of  values  when  it  is  present.  As  averaging  time  is  increased,  the  estimated 
range  rate  is  observed  to  coalesce  into  the  true  range  rate  (approximately  1.6  m/s  on  average),  even  when 
the  interferer  is  present.  This  is  accompanied  by  a  more  frequent  determination  of  the  correct  range  and 
depth  at  of  the  source  over  time. 

It  is  common  practice  to  average  ambiguity  surfaces  without  regard  to  source  motion  in  an  effort  to 
enhance  the  signal-induced  features  of  the  ambiguity  surface  relative  to  the  noise-induced  features.  Because 
of  the  broadness  of  the  peaks  in  Figs.  2  and  3  in  the  range-rate  direction  for  the  smaller  averaging  times, 
it  is  evident  that  some  track  detection  enhancement  may  be  realized  even  for  a  zero  range-rate  assumption. 
This  is  because  for  small  averaging  times  the  peaks  have  not  moved  enough,  relative  to  each  other,  to  cause 
significant  degradation  upon  averaging.  For  the  higher  averaging  times  (e.g.  128  samples),  however,  the 
consecutive  peaks  become  separated  enough  to  cause  degradation  upon  averaging.  In  Figs.  2  and  3,  this  is 
manifested  as  a  higher  resolution  in  range-rate,  so  that  the  zero-range  rate  line  no  longer  intersects  the  peaks. 
These  effects  are  illustrated  in  the  automated  results  in  Fig.  5,  which  plots  the  optimal  peak  data  versus 
times  for  a  zero  range  rate  (no  TBD)  assumption  for  64  and  128  samples.  It  is  observed  that  the  64-sample 
case  is  nearly  as  good  as  that  obtained  with  TBD  in  Fig.  4d.  This  is  because  the  source  peaks  have  not 
moved  much  relative  to  each  other  and  the  zero  range-rate  line  in  Figs.  2  and  3  intersect  the  maximum  peak. 
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For  128  samples,  however,  the  track  performance  without  TBD  in  Fig.  5b  has  been  destroyed,  because  the 
source  peaks  are  separated  enough  in  time  to  destroy  the  average  peak,  and  the  zero  range-rate  line  in  Figs. 
2  and  3  no  longer  intersect  the  maximum  peak.  When  track  motion  is  taken  in  account  using  this  averaging 
time  (Fig.  4e),  however,  the  good  tracking  performance  is  maintained.  The  maximum  correlation  versus 
time  plots  in  Fig.  4  include  a  dashed  curve  representing  the  value  obtained  without  TBD  (zero  range  rate). 
Note  the  significant  enhancement  in  correlation  for  the  128-sample  case. 

Little  is  known  firom  data  logs  or  other  sources  regarding  the  true  location  and  speed  of  the  interferer 
disrupting  the  signal  between  25  and  45  min  into  the  track.  It  is  interesting,  therefore,  to  let  the  tracker 
estimate  these  parameters  for  the  interferer.  This  was  done  by  applying  the  TBD  technique  to  data  in 
a  noise-only  bin  centered  at  96  Hz.  Figure  6  plots  the  automated  optimal  peak  results  versus  time  and 
averaging  time.  For  low  averaging  times,  it  is  diflScult  to  distinguish  any  clear  track;  the  data  is  largely 
scattered  over  all  values  of  range  rate,  range,  and  depth.  As  averaging  time  is  increased,  however,  the  range 
rate  is  seen  to  coalesce  toward  the  zero  range  rate,  except  during  the  period  of  the  interferer,  in  agreement 
with  the  behavior  observed  in  Figs.  2  and  3.  Range  and  depth,  however,  do  not  appear  to  follow  any 
preferential  pattern  when  the  interferer  is  absent,  suggesting  again  that  the  noise  during  these  periods  in 
probably  caused  by  long-range  shipping.  A  close  observation  of  the  128-sample  case  reveals  a  track  estimate 
for  the  interferer  has  been  successfully  obtained.  For  a  period  beginning  roughly  at  25  min  from  the  start  of 
the  track,  source  range  for  consecutive  times  steadily  increases  from  about  1  to  6  km,  source  depth  remains 
constant  at  about  30  m,  and  range-rate  congregates  around  10  m/s.  It  is  impossible  to  ascertain  the  accuracy 
of  these  estimates,  since  the  true  values  for  the  interferer  are  not  known,  and  since  the  interferer  was  most 
likely  in  a  range-dependent  environment  away  from  the  source  track.  In  addition,  the  peaks  tracked  may 
actually  be  sidelobes  of  the  true  interferer  peak.  Nevertheless,  these  results  demonstrate  the  feasibility  of 
this  approach  for  practical  scenarios. 

It  is  well  known  that  if  a  single  source  is  transmitting  a  signal  with  a  frequency  content  in  some  band 
of  frequencies,  then  ambiguity  surface  sidelobes  can  be  reduced  by  averaging  single-frequency  ambiguity 
surfaces  across  the  band  [5], [6], [7].  To  investigate  how  this  technique  effects  the  TBD  results,  the  shift- 
then-average  technique  was  applied  to  four-tone  (70,  95,  145,  and  195  Hz)  frequency-averaged  signal-fnoise 
surfaces  and  to  four-noise-oniy-bin  (71,  96,  146,  and  196  Hz)  frequency-averaged  noise-only  surfaces.  Results 
are  presented  in  Figs.  7  and  8,  in  the  same  format  as  Figs.  2  and  3.  Comparing  these  two  sets  of  figures, 
it  is  observed  that  the  range-rate  resolution  has  been  increased  for  the  frequency-averaged  surfaces,  since 
a  stricter  requirement  of  obtaining  a  match  for  all  four  frequencies  has  been  instituted.  In  addition,  it  is 
observed  that  the  track  is  now  much  more  detectable  during  the  period  when  the  interferer  is  present  for  all 
averaging  times.  This  is  particularly  striking  for  the  comparison  between  Fig.  3a  and  Fig.  8a.  This  track 
detection  enhancement  may  be  attributable  to  two  possible  influences.  First,  the  signal-to-interferer  ratios 
may  be  higher  for  the  additional  tonals,  resulting  in  less  interference  upon  averaging.  Second,  the  frequency¬ 
averaging  has  greatly  reduced  the  sidelobes  of  the  interferer,  allowing  the  source  to  be  more  detectable.  A 
comparison  of  the  noise-only  plots  in  Figs.  2  and  3  with  those  in  Figs.  7  and  8  reveals  that  the  responses  at 
high  range  rate  disappear  when  frequency  averaging  is  used,  suggesting  that  the  tracker  is  actually  tracking 
the  interferer’s  sidelobes  in  the  single-frequency  case.  The  observed  track  enhancement  is  most  likely  of 
combination  of  these  effects. 

Figure  9  presents  the  automated  optimal  track-parameter  estimates  for  the  frequency-averaged  case,  in 
the  same  format  as  Fig.  4.  It  is  observed  that  extremely  good  track  detection  and  parameter  estimation 
are  obtained  by  the  combination  of  TBD  and  frequency-averaging  techniques,  for  the  reasons  discussed 
above.  But  in  addition  to  a  reduction  of  the  interferer’s  sidelobes,  the  source’s  sidelobes  have  also  been 
reduced,  decreasing  the  possibility  that  a  source  sidelobe  will  be  mistaken  for  the  true  source  location. 
It  is  particularly  striking  that  very  good  track  estimates  can  be  obtained  with  small  averaging  times.  It 
should  be  noted,  however,  that  these  results  are  for  high  signal-to-noise  ratios,  at  least  when  the  interferer 
is  absent.  The  analysis  of  the  SWellEX-3  data  will  demonstrate  that  these  results  are  difficult  to  achieve  at 
low”  signal-to-noise  ratios. 
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IV.2  SWellEX«3  results 


The  SWellEX-3  data  analysed  here  was  characterized  by  a  low  signal-to~noise  ratio  over  the  entire  extent 
of  the  track.  During  the  first  30  min  of  the  80-min  data  set,  the  tow  ship  loitered  near  a  constant  position 
aw-aiting  commencement  of  the  source-tow.  As  the  source  moved  away  from  the  VLA  (with  a  speed  between 
2.4  and  2.9  m/s)  along  a  near-radial  track,  the  signal-to-noise  ratio  decreased  further,  providing  a  good  test 
on  the  limits  of  TBD  detection  performance.  Figures  10  and  11  present  maximum  correlation  in  averaged 
surfaces  versus  range  rate  and  time  for  this  data,  using  the  same  format  in  Figs.  2  and  3.  Results  for  three 
periods,  (a)  35  min,  (b)  45  min,  and  (c)  55  min  from  the  beginning  of  the  data  set,  are  included  in  Fig.  11 
so  that  the  performance  as  the  source  moves  away  from  the  VLA  can  be  observed.  Figure  10  indicates  that 
the  nearly  constant  position  of  the  source  during  the  loitering  period  is  clearly  detected.  However,  during 
the  period  of  the  source  tow,  the  track  is  diflScult  to  detect  at  the  lower  averaging  times.  The  track  at  a 
range-rate  of  about  2.5  m/s  (at  times  greater  than  30  min)  clearly  becomes  more  detectable  as  averaging 
time  is  increased,  and  is  accompanied  by  an  increase  in  range-rate  resolution.  Note  also  the  tendency  of 
noise-only  results  toward  a  zero  range  rate  over  the  entire  time  period,  suggesting  that  the  noise  in  this  data 
set  is  caused  primarily  by  long-range  shipping.  These  results  are  particularly  striking  in  Figs.  11a  and  11b. 
Note  in  both  Figs.  10  and  11,  however,  that  as  the  source  moves  further  away  from  the  VLA,  the  averaging 
times  considered  here  become  inadequate  for  track  detection.  For  example,  Figure  11c  fails  to  detect  a  track 
at  55  min  from  the  start  of  the  data  set  using  128  samples. 

Figure  12  presents  automated  track-parameter  estimates  for  optimal  tracks,  in  the  same  format  used  in 
Fig.  4.  As  with  the  SWellEX-1  data,  track  detection  and  parameter  estimation  improves  with  increasing 
averaging  time.  While  the  true  data  values  are  not  plotted  in  this  figure,  it  should  be  stated  that  the  detected 
tracks  agree  well  with  the  known  positions  of  the  source  during  this  period.  The  loitering  period,  extending 
out  to  about  30  min,  is  clearly  detected,  characterized  by  a  near-zero  range  rate,  a  nearly- constant  source 
range  (between  2  and  2.5  km),  and  a  nearly-constant  source  depth  of  about  50  m.  Beyond  the  loitering 
period,  the  2.5-m/s  source  tow  track  is  detected  out  to  a  time  of  about  50  min  when  128  samples  (5.87-min 
averaging  time)  are  used.  Note  that  the  source-tow  track  is  not  detected  using  less  than  32  samples.  Beyond 
50  min,  the  tracker  is  operating  on  the  ambient  noise,  with  range-rate  tending  toward  zero,  and  source  range 
and  depth  possessing  considerable  scatter. 

Results  for  the  case  in  which  TBD  techniques  are  not  implemented  and  the  surfaces  are  averaged  without 
searching  for  range  rate  (zero  range-rate  assumed)  are  presented  for  64  and  128  samples  in  Fig.  13.  Track 
detection  performance  is  seen  to  be  seriously  degraded,  failing  to  detect  the  source  tow  after  30  min.  The 
reason  this  failure  occurs  for  a  smaller  number  of  samples  than  observed  for  the  SWellEX-1  data  in  Fig. 
5  (64  instead  of  128)  is  the  fact  the  source  speed  was  larger  in  the  SWellEX-3  data  (2.5  m/s  instead  of 
1.6  m/s).  This  greater  speed  produces  a  greater  separation  between  consecutive  source  peaks,  resulting  in 
greater  degradation  upon  averaging  without  shifting.  Alternatively,  the  peak  in  the  maximum  correlation 
versus  range  rate  plots  of  Figs.  10  and  11  have  moved  further  to  the  right,  decreasing  the  tendency  of  the 
zero  range-rate  line  to  intersect  that  peak. 

The  effect  of  frequency  averaging  on  the  TBD  results  for  SWellEX-3  data  is  demonstrated  in  the  maximum 
correlation  versus  range  rate  curves  of  Fig.  14,  for  which  four  lovr  level  tonals  (77,  125,  157,  and  205  Hz) 
were  averaged.  The  curves  are  seen  to  be  sharper  and  cleaner,  as  a  result  of  the  added  requirement  that 
a  match  for  all  tonals  be  obtained.  The  lower  correlations  for  small  averaging  times  in  Figs.  14b  and  14c, 
compared  to  the  125-Hz  results  in  Fig.  11b  and  11c,  suggest  that  sidelobes  from  noise  sources  are  being 
reduced  by  the  frequency  averaging  process  at  these  times.  Nevertheless,  the  failure  to  detect  the  track  at 
55  min  from  the  start  of  the  data  set  (Fig.  14c)  implies  that  a  single  large  sidelobe  of  a  noise  source  is  not 
responsible  for  the  masking  the  of  source,  as  was  the  case  for  the  interferer  in  SWellEX-1.  Instead,  the  low 
level  of  the  source  has  caused  its  correlation  peak  to  lie  below  many  noise  sidelobes,  of  varying  level,  and 
frequency-averaging  is  unable  to  reduce  all  of  these  sidelobes  sufficiently  to  render  the  track  detectable. 

Finally,  the  automated  optimal  track-parameter  estimates  for  the  frequency-averaged  case  are  presented 
in  Fig.  15.  The  decreased  variability  of  the  estimates,  particularly  during  the  first  50  min  of  the  data  set,  is 
a  result  of  source  sidelobe  reduction,  which  decreases  the  possibility  that  a  source  sidelobe  will  be  mistaken 
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for  the  true  source  location.  But  once  again,  it  is  observed  that  frequency-averaging  fails  to  detect  the  track 
past  the  50-inin  limit. 


V  Detectibility 


A  measinre  of  the  matched-field  detection  performance  is  provided  by  the  detectability  (or  deflection  ratio), 
defined  as 


D  = 


10  log 


P signal+noise  l^noise 


(^noise 


(2) 


where  Psignai-h^noise  is  the  peak  level  in  a  range-depth  ambiguity  surface  (averaged  or  not)  containing  signal 
and  noise,  and  iinoUe  ^J^d  (Jnoise  the  mean  and  standard  deviation,  respectively,  of  a  range-depth  am¬ 
biguity  surface  containing  only  noise,  but  at  the  same  time  and  frequency  (or  close  to  that  frequency)  as 
the  signal+noise  surface.  Detectability  is  best  measured  experimentally  at  high  signal-to-noise  ratio,  when 
Psignai+noise  is  known  to  be  signal-induced.  When  noise  masks  the  signal  to  the  point  in  which  the  maximum 
peak  in  the  signal+noise  surface  is  noise-induced,  the  expression  above  no  longer  measures  the  detectability 
of  the  signal.  The  detectability  of  the  signal  in  low  signal-to-noise  ratio  environments  can  be  measured  if 
the  exact  source  location  is  known  and  Psignai^noise  above  is  replaced  by  the  correlation  at  that  location. 
In  this  paper,  however,  only  detectability  measurements  at  high  signal-to-noise  ratio  will  be  considered. 

A  high  signal-to-noise-ratio  case  suitable  for  calculating  detectability  is  that  following  the  period  of  the 
interferer  in  the  SWellEX-1  data.  Figure  16  plots  the  measured  detectability,  with  and  without  TBD,  versus 
the  number  of  samples  used  in  averages  beginning  52  min  from  the  start  of  the  track.  Also  shown  is  the 
UogM  theoretical  gain  expected  from  incoherently  averaging  the  surfaces;  on  the  log2  scale  used,  theory 
predicts  a  gain  of  1.5  dB  per  doubling  of  the  sample  size.  First  it  is  observed  that  when  source  motion  is  not 
taken  into  account  (zero  range-rate  assumption),  the  detectability  first  rises  gradually  according  to  theory 
and  then  drops  sharply  as  the  number  of  samples  is  increased.  The  drop  is  caused  by  the  source  peaks  for  a 
moving  source  not  being  aligned  when  averaging  over  time,  as  already  discussed.  When  shift-then-average 
TBD  technique  is  used,  however,  the  theoretical  gain  is  retained  because  the  TBD  technique  seeks  to  realign 
the  source  peaks  for  the  optimal  range-rate.  The  gain  in  detectability  for  the  TBD  case  is  about  8  dB 
above  that  obtained  without  TBD  when  128  samples  are  used.  The  higher-than-theoretical  gain  obtained 
using  TBD  at  the  larger  averaging  times  is  most  likely  caused  by  changes  in  source  or  noise  levels  over  the 
averaging  time. 


Conclusion 

The  feasibility  of  using  a  shift-then-average  TBD  technique  to  improve  detection  and  obtain  source  track 
parameters  has  been  experimentally  demonstrated  using  SWellEX-1  and  3  source-tow  data.  A  track  originally 
masked  by  a  strong  interferer  was  successfully  detected,  and  its  source  parameters  accurately  estimated, 
using  this  technique  in  the  SWellEX-1  data.  When  the  interferer  was  absent,  the  gains  in  detectability  were 
shown  to  agree  well  with  the  theoretical  value  expected  for  incoherently  averaging  the  surfaces.  Gains  in 
detectability  of  8  dB  over  a  that  obtained  by  averaging  without  regard  for  the  source  motion  were  realized. 
Previously  masked  tracks  were  also  detected  in  the  low  signal-to-noise  SWellEX-3  data  set,  which  was 
characterized  by  long-range  shipping  noise.  The  gain  with  TBD  over  that  without  TBD  was  greater  for  the 
SWellEX-3  data  because  of  the  greater  speed  of  the  source-tow.  Frequency  averaging  appears  to  improve 
TBD  performance  primarily  in  cases  in  which  the  source  peaks  are  masked  by  strong  sidelobes  of  noise 
sources,  as  occurs  when  a  strong  interferer  passes  nearby. 

An  extension  of  this  technique  to  the  problem  of  detecting  short-range  off-radial  tracks  is  currently 
underway. 
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Figure  4.  TBD  (shift-then-average)  automated  detection  results  for  95-Hz  tonal  in  SWeUEX-l 
source-tow  track.  Maximum  correlation  and  range,  depth,  and  speed  (range  rate)  at  maximum 
are  plotted  versus  time  for  (a)  8  samples  (0.32  min),  (b)  16  samples  (0.68  min),  (c)  32  samples 
(1.41  min),  (d)  64  samples  (2.87  min),  and  (e)  128  samples  (5.78  min)  in  average.  The  dashed 
curves  in  the  maximtim  correlation  versus  time  plots  are  those  obtained  without  TBD,  as  in 
Fig.  5. 
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Figure  5.  Automated  detection  results  obtained  without  TBD  (averaging  ambiguity  surfaces 
without  searching  for  range  rate)  for  the  95-Hz  tonal  in  SWellEX-1  source-tow  track.  Maximum 
correlation,  range  at  maximum  correlation,  and  depth  at  maximum  correlation  are  plotted  versus 
time  for  (a)  64  samples  and  (b)  128  samples  in  average.  The  speed  (or  range  rate)  in  the  last 
column  is  always  zero  for  this  case. 
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Figure  6.  TBD  (shift-then-average)  automated  detection  results  for  96-Hz  noise-only  bin  in 
SWellEX-1  source-tow  track.  Maximum  correlation  and  range,  dep±,  and  speed  (range  rate) 
at  maximum  are  plotted  versus  time  for  (a)  8  samples  (0.32  min),  (b)  16  samples  (0.68  min), 

(c)  32  samples  (1,41  min),  (d)  64  samples  (2.87  min),  and  (e)  128  samples  (5.78  min)  in  average. 
The  dashed  curves  in  the  maximum  correlation  versus  time  plots  are  those  obtained  without  TBD. 
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Hgure  7.  Maximum  correlation  in  frequency-averaged  shifted-then-averaged  range-depth 
ambiguity  surfaces  versus  candidate  range  rate,  time,  and  number  of  samples  in  average  for 
SWellEX-1  source-tow  track.  Number  of  samples  in  average  =  (a)  16  (0.68  min),  (b)  32 
(1.41  min),  (c)  64  (2.87  min),  and  (d)  128  (5.78  min).  First  column  is  for  a  four-tone  (70, 95, 
145,  and  195  Hz)  frequency  average,  while  second  column  is  for  a  four  noise-only-bin  (71, 
96, 146,  and  1%  Hz)  frequency  average. 
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Figure  9,  TBD  (shift-then-average)  automated  detection  results  for  four-tone  frequency  average 
in  SWelIEX-1  source-tow  track.  Maximum  correlation  and  range,  depth,  and  speed  (range  rate) 
at  maximum  are  plotted  versus  time  for  (a)  8  samples  (0.32  nain),  (b)  16  samples  (0.68  min),  (c) 
32  samples  (1.41  min),  (d)  64  samples  (2.87  min),  and  (e)  128  samples  (5.78  min)  in  average. 

The  dashed  curves  in  the  maximum  coixelation  versus  time  plots  are  those  obtained  without  TBD. 
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Figure  1 L  Maximum  correlation  in  shifted-then-averaged  range-depth  ambiguity  surfaces 
versus  candidate  range  rate  and  number  of  samples  in  average  at  three  periods  of  SWellEX-3 
source-tow  track.  Averaging  begins  at  (a)  35  min,  (b)  45  min,  and  (c)  55  min,  from  the  start 
of  the  track,.  In  each  row,  the  first  column  is  for  the  125-Hz  tonal  (S+N)  while  the  second 
column  is  for  a  noise-only  bin  centered  at  126  Hz.  The  vertical  dashed  lines  mark  the  true 
range  rate  of  the  source-tow  during  each  period. 
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Figure  14.  Maximum  correlation  in  frequency-averaged  shifted-then-averaged  range-depth 
ambiguity  surfaces  versus  candidate  range  rate  and  number  of  samples  in  average  at  three 
periods  of  SWellEX-3  source-tow  track.  Averaging  begins  at  (a)  35  min,  (b)  45  min,  and  (c) 
55  min,  from  start  of  track.  Frequency  average  is  over  four  tones  (77, 125, 157,  and  205  Hz). 
The  vertical  dashed  lines  mark  the  true  range  rate  of  the  source-tow  during  each  period. 
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Figure  15.  TBD  (shift-then-average)  automated  detection  results  for  four-tone  (77, 125, 157, 
and  205  Hz)  frequency  average  in  S WellEX-3  source-tow  track.  Maximum  correlation  and 
range,  depth,  and  speed  (range  rate)  at  maximum  are  plotted  versus  time  for  (a)  8  samples 
(0.32  min),  (b)  16  samples  (0.68  min),  (c)  32  samples  (1.41  min),  (d)  64  samples  (2.87  min), 
and  (e)  128  samples  (5.78  min)  in  average.  The  dashed  curves  in  the  maximum  correlation 
versus  time  plots  are  those  obtained  without  TBD. 


No.  of  Samples 

Figure  16.  Detectibility  versus  number  of  samples  (averaging  time)  for  averages  beginning 
at  52  min  from  start  of  track  in  SWelIEX-1 .  Solid  line:  TBD  (shift-then-average)  result.  Top 
dashed  line:  5  log  M  theoretical  gain  curve,  where  M  =  number  of  samples.  Bottom  dashed  line: 
No  TBD  (averaging  without  range-rate  search)  result. 


51 


Matched-Field  Processing  in  a  Range-Dependent 
Shallow  Water  Environment  in  the  North-East  Ocean: 

Array  Tilt  Considerations 

M.  L.  Yeremy,  J.  M.  Ozard, 

Defence  Research  Establishment  Atlantic 
Esquimalt  Defence  Research  Detachment 
FMO  Victoria,  B.C.,  Canada  VOS  IBO 

N.R.  Chapman 
University  of  Victoria 
Victoria,  B.C.,  Canada,  V8W  2Y2 

and 

M.J.  Wilmut 
Royal  Military  College 
Kingston,  Ontario,  Canada  K7K  5L0 

3  SEPTEMBER  1996 


Introduction 

A  series  of  ocean  acoustic  experiments  referred  to  as  the  PACIFIC  SHELF  Sea  Trial  was  completed  in 
September,  1993  by  the  Defence  Research  Establishment  Pacific  (DREP),  Victoria,  B.C.  and  the  Applied 
Research  Laboratory  (ARL),  University  of  Texas  at  Austin.  The  purpose  of  these  experiments  was  primarily 
to  evaluate  Matched-Field  Processing  (MFP)  techniques  in  a  Pacific  shallow  water  environment.  Efficient 
Linear  Tracker  (ELT)  results  firom  this  trial  are  shown  in  this  paper,  for  data  collected  from  a  Vertical  Line 
Array  (VLA)  during  an  experiment  where  a  Continuous  Wave  (CW),  multi-frequency  source  was  towed  from 
a  shallow-water  position  on  the  shelf  towards  the  VLA,  farther  downslope.  In  particular,  the  effect  that  array 
tilt  has  on  these  results  is  investigated. 


I  Environment,  Data  and  Processing 

The  towed  source’s  track,  shown  in  Figure  1,  began  11.9  km  firom  the  array,  on  the  continental  shelf  where 
the  water  column  depth  was  ~  150  m  and  proceeded  down  the  continental  slope  towards  the  VLA  at  a 
source  depth  of  ~  30  m.  The  VLA  consisted  of  16  phones  which  were  equispaced  at  15  m  and  spanned  the 
depths  from  90  m  to  315  m  at  a  water  column  depth  of  ~  375  m.  Both  the  source  and  receiver  positions 
were  recorded  from  GPS  (Global  Positioning  System)  and  Radar  which  each  have  range  errors  of  ~  200  m. 
The  experiment  was  conducted  during  a  time  when  there  was  considerable  shipping  traffic. 

The  multi-frequency,  CW  source  emitted  three  tones,  45,  70  and  72  Hz,  each  of  a  different  and  specified 
Source  Level  (SL).  The  source  levels  at  45  Hz  and  70  Hz  were  typical  of  a  strong  line  on  a  merchant  vessel 
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Figure  1:  The  line  overlaying  the  region’s  bathymetry,  represents  the  towed  source’s  track  with  the  tow 
direction  shown  by  the  arrow.  The  star  shows  where  the  VLA  was  located  at  the  beginning  of  the  experiment. 


Sound  Speed  (m/s) 


Figure  2:  The  sound-speed  profiles  used  in  the  environmental  model  are  shown  as  well  as  the  hydrophone 
depths.  The  shear  speed  (shaded  and  dashed)  and  the  compressional  speed  (dark  and  solid)  axe  shown  on 
the  two  lower  abscissa  scales. 
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while  the  72  Hz  line  was  20  dB  lower.  The  results  obtained  for  a  strong  (45  Hz)  and  a  weak  (72  Hz)  line  are 
presented  here. 

The  environment  model  used  here  is  shown  in  Figure  2.  It  was  based  on  an  XSV  (Expendable  Sound 
Velocimeter)  measurement  which  was  cast  eighty  minutes  prior  to  the  experiment’s  start  time,  and  typical 
bottom  acoustic  values  (e.g.  Ref.  [1]  and  [2])  for  the  vicinity  of  the  array. 

Recent,  seismo^acoustic  results  from  another  PACIFIC  SHELF  experiment  [3]  indicate  that  the  sediment 
thickness  and  compressional  speed  of  the  bottom  half-space  in  the  location  of  the  VLA  are  similar  to  the 
values  used  for  this  analysis.  It  is  likely  that  these  bottom  parameters  vary  from  the  continental  shelf  to  the 
slope  regions.  However,  in  this  study,  we  did  not  have  suflBcient  knowledge  of  the  environment  to  model  the 
range-dependent  acoustic  parameters.  Therefore,  with  the  exception  of  the  bathymetry  shown  in  Figure  1, 
the  environment  model  used  here  was  the  same  at  all  ranges.  More  detailed  information  of  the  environment 
model  and  the  data  axe  given  in  [4]. 

The  results  presented  here  are  from  an  ELT  [5]  which  finds  constant  depth  linear  tracks  and  estimates 
their  track  Signal-to-Noise  Ratios  (SNR)  from  a  collection  of  MFP  ambiguity  surfaces  which  are  sampled  in 
time.  The  ELT  operated  on  a  Bartlett  ambiguity  surface  and  is  described  further  in  [4]. 

The  adiabatic  normal  mode  approximation  was  the  method  used  to  calculate  the  range-dependent  replica 
fields  required  for  MFP.  One  advantage  of  using  this  method  is  that  field  calculations  are  fast  and  efficient  if 
the  acoustic  modal  data  are  precomputed.  This  approach  does  not  include  effects  due  to  mode-coupling,  and 
inaccuracies  commonly  associated  with  these  effects  for  a  range-dependent  environment  were  addressed  in 
[4]  for  this  experiment.  It  was  found  that  mode-coupling  inaccuracies  were  negligible  except  at  the  steepest 
slopes  of  the  environment.  Nevertheless,  the  MFP  results  had  very  high  correlations  with  average  values 
of  about  0.85  and  0.70  for  the  respective  45  and  72  Hz  data  and  with  correlations  as  high  as  ^  0.95  for 
the  stronger  received  signals  at  45  and  70  Hz.  This  indicates  that  mismatch  from  the  replica  field  does 
not  substantially  debilitate  the  MFP  results  despite  the  mode-coupling  inaccuracies  introduced  with  the 
adiabatic  normal-mode  method. 

ORCA  [6],  [7],  [8]  was  the  normal  mode  program  used  for  modelling  the  environment’s  modal  data.  It 
was  used  primarily  because  the  data  were  collected  in  a  shallow  water  environment  which  hkely  had  low  sea 
floor  shear  speed  properties. 

In  Figure  1,  it  can  be  seen  that  the  ship  track  was  not  along  a  radial  path  relative  to  the  VLA  but 
consisted  of  two  segments  with  path  directions  on  different  headings.  For  this  study,  the  replica  field  was 
approximated  as  a  radial  instead  of  modelling  and  tracking  in  3-dimensions.  Because  a  linear  tracker  assumes 
a  constant  velocity  source,  the  data  were  analysed  separately  for  the  two  segments  which  are  distinguished  by 
the  heading  of  the  towed  source.  The  fax  and  near-range  segments  refer  to  the  furthest  and  nearest  segment 
relative  to  the  VLA. 


II  Results 

11,1  Performance  Assuming  No  Array  Tilt 

The  two  segments  were  processed  with  the  ELT  algorithm  assuming  source  depths  between  10  and  100  m 
in  10  m-depth  increments.  The  maximum  track  SNR  at  both  frequencies  for  both  segments  is  plotted  in 
Figures  3  and  4  respectively.  The  SNR  threshold  for  a  Probability  of  False  Alarm  of  ~  10“^  at  any  specific 
depth  is  found  to  have  a  linear-scale  SNR  level  of  8  (i.e.,  9  dB)  when  the  noise  is  spatially  white  [9],  [5]. 
A  signal  of  9  dB  SNR  would  have  a  probabihty  of  detection  of  0.75.  In  this  analysis,  the  source  track  is 
identified  by  the  largest  track  SNR  over  all  depths  which  is  also  greater  than  this  threshold.  The  source 
track  SNR  levels  for  the  72  and  45  Hz  far-range  segment  were  respectively  34  and  107  (15.3  and  20.3  dB), 
and  for  the  near-range  segment  these  were  98  and  661  (19.9  and  28.2  dB).  The  information  presented  in 
Figures  3  and  4  indicates  that  the  tracks  from  both  frequencies  are  well  above  the  threshold  SNR  level  and 
that  they  accurately  identify  the  source  depth  except  for  the  45  Hz  fax-range  segment  for  which  a  depth 
10  m  too  shallow  was  indicated.  For  this  array  design,  spatial  leakage  causes  the  track  SNR  to  be  above  the 
threshold  at  non-source  depths.  However  the  maximum  SNR  occurs  within  10  m  of  the  known  source  depth. 
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Figure  5:  The  45  and  72  Hz  track  results  are  shown  in  the  respective  left  and  right-hand  side  plots.  The 
dashed  and  dotted  lines  are  the  top  acoustic  tracks  from  the  respective  far  and  near-range  segments  and  the 
solid  line  is  the  GPS  track. 
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Figure  6:  Maximum  track  SNR  at  45  Hz,  for  the  far  (solid  hne)  and  near-range  (dashed  hne)  segments  as  a 
function  of  array  tilt. 


Figure  7:  Maximum  track  SNR  at  72  Hz,  for  the  far  (solid  line)  and  near-range  (dashed  line)  segments  as  a 
function  of  array  tilt. 


Figure  5  shows  the  range-time  tracks  for  the  45  Hz  and  72  Hz  signals.  The  most  significant  track  estimates 
for  both  data  segments  (at  the  depth  of  the  maximum  SNR  in  Figures  3  and  4),  are  plotted  as  dashed  and 
dotted  hnes  while  the  target  track,  which  is  based  on  GPS  measurements,  is  represented  as  a  solid  line. 
The  most  significant  tracks  for  both  the  45  and  72  Hz  signal  are  essentially  coincident  with  the  GPS  track. 
Variations  between  the  GPS  and  acoustic  tracking  are  at  most  300  m  and  the  maximum  difference  between 
the  72  and  45  Hz  estimated  tracks  was  130  m.  The  agreement  between  GPS  and  acoustic  tracks  at  both 
frequencies  is  remarkably  good  considering  that  many  other  sources  were  present  in  the  region.  In  this 
respect,  it  should  also  be  noted  that  the  SNR  at  72  Hz  is  low. 

It  should  be  noted  that  at  the  depths  with  the  largest  track  SNRs  the  top  25  track-positions  and  SNR 
are  very  similar  to  each  other.  At  depths  with  lower  track  SNR  the  most  significant  track  positions  are  much 
more  ambiguous  since  noise  and  depth  mismatch  affect  these  weaker  sidelobe  tracks  to  a  greater  extent. 

11,2  Performance  for  Diflferent  Array  Tilts 

A  MFP  mismatch  study  [10],  [11]  found  that  array  tilt  was  a  significant  source  of  mismatch.  In  this  study, 
tracking  was  performed  for  different  array  tilts  in  the  plane  of  the  radial.  The  results  are  shown  in  Figures  6 
and  7  for  the  respective  45  and  72  Hz  data.  The  track  SNR  is  consistent  with  an  array  tilt  of  ^  —2.5°  ±  0.5° 
as  the  SNR  improved  by  as  much  as  ~  2.5  dB  for  the  72  Hz  data  and  1.0  dB  for  the  45  Hz  data  compared 
to  no  array  tilt.  Here,  a  positive  tilt  implies  that  the  bottom  array  element  is  closest  to  the  source.  The  72 
Hz  track  estimates  were  nearer  to  the  GPS  ranges  by  as  much  as  120  m  while  the  45  Hz  results  improved  by 
30  m  over  the  whole  track  for  a  —2.5°  tilt.  Higher  correlations  were  also  observed.  Similar  array  tilts  (2.6°) 
were  measured  near  the  time  of  the  experiment.  However,  the  orientation  of  the  array  tilt  was  unknown. 


Ill  Conclusions 

MFP  with  tracking  was  applied  to  data  from  a  north-east  Pacific  shallow  water  environment.  Despite 
incomplete  environmental  knowledge  of  the  geoacoustic  bottom  profile,  high  MFP  correlations  were  obtained. 
The  average  of  the  largest  correlations  for  high  SNR  data  at  45  Hz  and  low  SNR  data  at  72  Hz  were  about 
0.85  and  0.70  respectively.  These  results  are  very  encouraging  since  they  imply  that  detailed  environmental 
knowledge  may  not  be  a  prerequisite  for  source  localization  with  MFP.  Better  environmental  information 
would  likely  enable  higher  track  SNRs  and  result  in  improved  source  tracking.  This  has  been  demonstrated 
here,  to  an  extent,  for  the  case  where  array  tilt  considerations  improved  the  SNR. 
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An  Efficient  Linear  Tracker  which  finds  the  most  significant  tracks  (over  the  set  of  ambiguity  surfaces) 
was  found  to  agree  with  the  GPS  source  range  within  300  m  and  known  source  depth  within  10  m  at  both 
frequencies.  The  accuracy  of  these  track  ranges  was  comparable  to  the  GPS  source  range  error  of  200  m. 
The  track  SNRs  for  the  72  and  45  Hz  far-range  segment  were  respectively  12,7  and  19.2  dB  while  for  the 
near-range  segment  these  were  respectively  17.7  and  28  dB,  The  SNR  and  track  estimates  were  improved  by 
as  much  as  2.5  dB  and  120  m  respectively,  when  array  tilt  was  taken  into  account.  The  highest  SNR  levels 
attained  occurred  for  a  modelled  array  tilt  which  was  very  near  recorded  array  tilt  values.  When  this  tilt 
was  modelled,  the  estimated  track  was  even  closer  to  the  GPS  track.  In  all  cases  the  tracker  contributed  to 
system  gain  due  to  the  integration  of  MFP  ambiguity  surfaces  over  time,  despite  the  presence  of  mismatch. 
These  encouraging  results  may  be  improved  by  a  better  knowledge  of  the  environment,  3D  modelling  and 
the  use  of  multiple  frequencies  in  the  analysis. 
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Abstract 

An  efficient  technique  for  determining  the  performance  of  Matched  Field  Processing  (MFP)  with 
tracking  is  summarized  and  illustrated  with  examples.  In  this  algorithm  the  candidate  source  tracks 
Me  the  linear  tracks  passing  through  the  strongest  peaks  in  ambiguity  surfaces  from  MFP  at  different 
times.  The  candidate  track  with  the  largest  Bartlett  sum,  over  the  sequence  of  ambiguity  surfaces  in 
time,  is  declared  a  source  track,  provided  it  is  above  a  preassigned  threshold.  Simulation  and  an  analytic 
approximation  are  used  to  determine  the  probability  the  source  track  is  examined.  Furthermore,  the 
probability  the  source  track  is  detected,  if  examined,  is  given.  The  above  probabilities  are  used  to 
illustrate  the  tracking  performance  of  slanted  and  horizontal  arrays  operating  in  an  Arctic  environment. 


Linear  Tracking  Algorithm 

Matched  Field  Processing  (MFP)  is  an  advanced  signal  processing  method  for  the  localization  and  detection 
of  acoustic  sources.^  In  this  paper  the  problem  of  detection  of  sources  of  low  signal- to-noise  ratio  (SNR)  is 
considered.  Signals  from  such  sources  are  matched  against  predictions  of  the  received  signal  for  ail  possible 
positions  to  form  ambiguity  ‘surfaces’.  A  Bartlett  beamformer  was  chosen  for  the  matching  in  this  study. 
The  SNR  is  so  low  that  the  Bartlett  statistic  at  the  source  position  is  often  not  the  largest  value.  Tracking, 
combining  the  information  on  a  set  of  contiguous  surfaces  is  used  to  determine  the  source’s  track,  if  one  is 
present.  In  this  study  we  assumed  a  source,  if  present,  is  moving  linearly  at  constant  speed,  depth  and  head¬ 
ing  and  that  NS  contiguous  ambiguity  surfaces  are  available  for  analysis.  The  unweighted  Linear  Tracking 
Algorithm  (LTA)  consists  of  five  sets  of  computations  performed  for  each  possible  depth: 

(1)  for  each  of  the  NS  surfaces  the  positions  of  the  NPK  largest  peaks  are  determined; 

(2)  all  possible  linear  tracks  joining  pairs  of  peaks  on  different  surfaces  are  found.  These  are  called  combi¬ 
natoric  tracks; 

(3)  a  constraint  to  realistic  maximum  speeds  for  the  source  is  imposed  to  reduce  the  combinatoric  tracks  to 
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the  physically  possible  tracks; 

(4)  the  unweighted  track  statistic  is  calculated  for  each  physically  possible  track: 

.  NS 

where  B{Ti)  is  the  Bartlett  output  for  the  unit  norm  replica  vector  r,  on  the  surface  along  the  track; 

(5)  the  unweighted  estimated  track  SNR  outputs  are  calculated: 

(2) 

S 

where  x  is  the  average  value  and  s  is  the  standard  deviation  of  the  Barlettt  statistic  for  all  points  on  all 
ambiguity  surfaces  in  noise  alone. 

To  be  detected  a  source  track  must  be  examined  (i.e.  be  one  of  the  tracks  found  at  step  2),  and  if  examined 
have  the  largest  estimated  track  SNR  and  exceed  a  threshold.  As  the  Bartlett  statistics  on  a  surface  are 
correlated  the  probabilities  for  the  above  two  events  are  difficult  to  determine.  A  technique  employing  some 
simulation  and  an  analytic  approximation  will  be  described  below  to  determine  the  Probability  the  source 
Track  is  Examined  (PTE).  Reference[5]  contains  results  needed  to  calculate  the  Probability  the  source  Track 
is  Detected,  if  Examined  (PTDE). 


Simulation  Scenario 

A  range-independent  upward-refracting  channel  with  a  650  m  water  depth  and  a  22-m  thick  attenuating 
sediment  layer,  representative  of  an  Arctic  scenario,  was  chosen.^  A  normal  mode  propagation  model  was 
used  to  model  the  25  Hz  sound  source.  The  sound  source  was  at  a  depth  of  100  m  and  spatially  white  noise 
was  added  to  the  signal.  Data  was  simulated  for  a  20  element  slanted  array  (vertically  spacing  spanning  the 
water  column  and  with  half  wavelength  horizontal  spacing),  a  10  element  slanted  array  (spanning  the  top 
half  of  the  water  column  and  with  half  wavelength  horizontal  spacing),  and  a  20  element  horizontal  array 
(located  on  the  bottom  with  half  wavelength  spacing).  The  arrays  were  oriented  so  their  horizontal  extent 
was  along  the  x  axis.  The  search  region  was  two  dimensional,  5  to  30  km  in  range  and  -45®  to  45®  in  bearing. 
K  averages  were  used  in  the  formation  of  the  data  covariance  matrix  used  in  the  Bartlett  processing. 


Probability  Source  Track  Examined 

The  probability  the  source  track  is  examined  is 

NS  NS  NS 

prE  =  i-[I][(i-P.)  +  E(P‘)  n  (1-Pi)]  (2) 

»=1  *=1 


where  pi  is  the  probability  that  the  Bartlett  statistic  for  a  region  about  the  source  forms  one  of  the  top 
NPK  peaks  on  the  surface.  This  region  must  be  sufficiently  large  to  include  source  peaks  that  have  been 
displaced  by  noise  effects  and  yet  not  too  large  so  as  to  include  neighbouring  non  source  peaks.The  pi  are 
dependent  on  many  factors  and  can  be  found  by  simulation.  In  Reference  4  it  is  shown  that  if  the  region 
about  the  source  is  0.15-0.20  of  the  signal  half  width  (a  source  position  dependent  quantity)  then  pi  may  be 
approximated  by  PT,R{sy 


Pt,R{9)  ” 


1 


r 


e  ^  dt 


(4) 


Pt,R{s)  is  the  probability  that  the  Bartlett  output  for  source  position  vector  s  is  greater  than  the  threshold 
T.  iZ(s)  is  the  source  SNR  in  the  ambiguity  surface.  T  is  the  threshold  above  which  there  will  be  NPK  peaks 
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in  the  noise-only  standarized  Bartlett  ambiguity  surface.  The  threshold  is  a  function  of  the  scenario,  the 
number  of  peaks  ranked  and  the  number  of  averages  K.  The  per  sensor  SNR  is  assumed  to  be  small  and  K 
is  assumed  to  be  large.  The  threshold  values  for  NPK  peaks  are  found  by  simulation.  For  any  track  PTE 
was  found  from  Equation  3  using  the  approximation  to  the  pi’s  along  the  track  from  Equation  4. 


Array  comparison 

The  performance  of  the  above  arrays  was  found  for  a  source  sequentially  positioned  at  each  of  five  equispaced 
locations  along  a  linear  track.  The  initial  position  was  10000  m  range  at  a  -30^  bearing  and  the  final  position 
was  25000  m  range  at  a  30®  bearing-  The  source  level  was  105  dB  re  1  fiPa  at  1  m  and  the  noise  level  was 
50  dB  re  1  ^Pa  at  all  sensors.  Figure  1  shows  the  threshold  T  versus  the  number  of  peaks  NPK  and  the 
number  of  averages  K  for  the  three  arrays.  PTE  and  PTDE  are  tabulated  in  Table  I  for  the  three  arrays  for 
an  NPK  of  20  and  a  K  of  60. 
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Figure  1:  Estimated  thresholds  for  peaks  in  the  noise-only  standarized  ambiguity  surface  as  a  function  of 
number  of  peaks  ranked,  for  the  number  of  averages  indicated:  (a)  slanted  20-element  array;  (b)  slanted 
10-element  array;  (c)  horizontal  20-element  array.  The  threshold  scale  is  linear  in  terms  of  the  standard 
deviation  of  the  noise. 

Not  suprisingly  PTE  is  high  when  PTDE  is  high  and  visa-versa.  As  can  be  seen  the  20  element  slanted 
array  would  very  likely  detect  the  track,  the  10  element  slanted  array  might  detect  the  source  track  while 
the  horizontal  array  would  almost  never  detect  the  source.  The  difference  in  performnce  between  the  arrays 
stems  from  the  difference  in  SNR  received  at  the  arrays  through  the  term  i?(s)  in  Equation  4.  This  study 
does  not  adjust  the  noise  level,  which  would  probably  be  significantly  lower  for  the  bottom  horizontal  array, 
and  thus  dramatically  underestimates  the  horizontal  arrays  performance.  This  could  be  taken  into  account 
in  future  work. 
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Table  1:  Tracking  performance  of  the  three  arrays  for  the  above  scenario.  The  probability  the  track  was 
examined  (PTE)  and  the  probability  the  source  track  is  detected,  if  examined  (PTDE).  The  threshold  for 
PTDE  corresponds  to  a  false  alarm  probability  of 


Array  Type 

PTE 

PTDE 

Slanted  20  element 

0.974 

0.90 

Slanted.  10  element 

0.892 

0.60 

Horizontal  20  element 

0.005 

0.01 

Concluding  Remarks 

The  probability  a  source  track  is  examined  in  a  particular  scenario  using  the  Linear  Tracking  Algorithm  may 
be  found  with  minimal  computer  simulation  and  an  analytic  approximation.  The  probability  a  source  track, 
if  examined,  is  detected  may  also  be  determined  analytically.  This  approach  may  be  applied  to  any  array 
configuration  to  estimate  array  detection  performance. 
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Broad-band  Shallow  Water  Detection  Experiments: 
Adaptive  Matched-Field  Processing  on  VLA's  and 
Adaptive  Plane-Wave  Beamforming  on  HLA's 

Newell  Booth  and  Phil  Schey,  NRaD 
William  Hodgkiss,  SIO/MPL 


Abstract 

Quantitative  measurements  of  multi-tone  detection  performance  are  presented  for  two 
broadband  adaptively  processed  array  systems: 

1)  vertical  line  arrays  with  matched  field  rephcas,  and 

2)  horizontal  line  arrays  with  plane  wave  beamforming  replicas. 

The  issues  of  measuring  minimum  detectable  level,  array  gain,  and  deflection  ratios  are  discussed 
and  compared.  Clutter  rejection  and  target  classification  by  depth  are  illustrated  for  the  matched 
field  processed  VLA.  Remaining  critical  issues  for  use  of  matched  field  processing  in  undersea 
surveillance  are  summarized. 

The  results  presented  were  obtained  in  the  third  Shallow  Water  Evaluation  Cell 
Experiment  (SWellEX-3)  was  conducted  12  km  off  the  coast  of  San  Diego,  CA,  in  July  1994'-^. 
The  conditions  of  the  experiment  are  shown  in  Fig.  1 .  Data  were  recorded  on  a  vertical  line  array 
(VLA)  suspended  from  the  research  platform  FLIP  and  on  a  bottom  mounted  horizontal  line  array 
(HLA)  at  the  positions  shown. 

SWelIEX-3  Experiment,  ArcMFP-5  Event 

On  30  July  between  13:10  and  13:30  GMT,  an  acoustic  source  (J-13)  was  towed  at  30  m 
depth  along  the  track  shown  in  Fig.  1.  Propagation  from  the  source  was  along  the  200  m  isobath. 
The  sound  speed  profile  measured  at  FLIP  is  also  shown.  The  acoustic  source,  towed  at  30  m 
depth,  transmitted  a  continuous  50-tone  comb  signal,  illustrated  on  the  right  side  of  Fig.  1. 

This  multi-tone  Minimum  Detectable  Level  (MDL)  signal  consisted  of  ten  pilot  tones 
spread  between  53  and  197  Hz.  Each  pilot  tone  at  frequency  f  -,  was  accompanied  be  a  set  of  four 
tones  at  frequencies  y^+2,/+4,yj+6,  and yj+8  Hz,  at  levels  -10,-14,-18  and  -22  dB  below  the  pilot 
tone,  respectively.  The  results  reported  in  this  paper  use  the  pilot  tones  to  estimate  the  input 
signal  to  noise  ratio  (ISNR)  from  which  the  ISNR  of  the  lower  level  tones  is  inferred  by 
subtracting  the  appropriate  difference  in  source  level.  The  ISNRs  for  the  pilot  tones  and  the  “-22 
dB”  set  of  tones  is  shown  in  the  lower  left  of  Fig.  1  for  event  ArcMFP-5.  This  signal  enables 
calibrated  measurement  of  the  performance  of  the  adaptive  algorithms  at  low  signal  to  noise  ratio, 
where  the  algorithms  are  less  sensitive  to  replica  mismatch. 
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Beamforming  and  Adaptive  Processing  Methods 

Processing  was  done  using  the  ten  pilot  tone  frequencies,  the  “-22  dB”  set  and  a  noise  I 

only  set  of  ten  frequencies  at  fj+9  Hz,  1  Hz  above  the  “-22  dB”  set.  Adaptive  Matched-Field 
Processing  (MFP)  was  done  on  a  78.5  m ,  =8  element  subset  of  the  64  available  array  elements  I 

of  the  VLA  QJl  @  67  Hz).  Adaptive  Plane- Wave  Beamforming  (PWB)  was  also  done  on  a  52  I 

m,  =14  element  subset  of  the  HLA  (7/2  @  187.5  Hz).  Each  six  seconds  of  element  data  was  J 

Fourier  analyzed  with  a  50%  overlapped  Kaiser-Bessel  window  to  form  a  element  data  vector,  I 

X  ,  for  each  pilot  frequency,/  ,  and  time,  .  The  frequencies  of  the  pilot  tones  were  ■ 

tracked  and  used  to  derive  the  low  level  and  noise  only  frequencies  that  were  used  in  the 
processing.  The  covariance  matrix  was  estimated  by  averaging  N,  outer  products  of  the  data  I 

vectors  * 

1  i 

KW.O  =  /  Ex(y;..(,-)X‘(y;.,r)  (1)  I 

A,  y=i  I 

I 

I 
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for  T  seconds,  where  *  indicates  complex  conjugate  transpose.  The  processing  and  array 
parameters  are  summarized  in  TABLES  I  and  II  for  the  VLA  and  HLA  respectively. 

TABLE  I:  VERTICAL  ARRAY  AND  PROCESSING  PARAMETERS 


Name 

Value 

Comments 

Number  of  elements, 

8 

Spacing  1 1.25  m  (A/2  at  67  Hz) 

Array  Aperture 

78.5  m 

starting  6  m  firom  bottom  at  198  m  depth 

Array  Tilt  (a) 

<2  deg 

Toward  North  during  both  events 

Sampling  rate 

1,500  Hz 

FFT  length 

8192 

50  %  overlapped 

FFT  binwidth 

0.18  Hz 

Kaiser-Bessel  weighting  with  a  =  2.5 

Covariance  Integration  Time,  T 

60  sec 

Nt  =  20  samples 

TABLE  II:  HORIZONTAL  Array  and  Processing  Parameters 


Name 

Value 

Comments 

Number  of  elements, 

14 

Spacing  4  m  (A/2  at  187.5  Hz) 

Array  Aperture 

52  m 

Sampling  rate 

1,562.5  Hz 

FFT  length 

8192 

50  %  overlapped 

FFT  binwidth 

0.19  Hz 

Kaiser-Bessel  weighting  with  a  =  2.5 

Covariance  Integration  Time,  T 

73  sec 

N,  =  28  samples 

The  narrowband  beamforming  (matched- field  or  plane- wave)  process  used  was  an 
estimate  of  the  normalized  cross-correlation  of  the  received  data  with  weight  vectors,  w.  For 
each  time,  t,  and  pilot  fi^equency,/,  the  “correlation”  was  calculated  using 


w-(y;.,r)  K{f.,t)  w(y;.,r) 


(2) 


where  Tr(K)  is  the  trace  of  the  covariance  matrix,  K.  Since  Tr(K(X.,»))  =  iV/A»W).<he 
correlation  is  the  output  power  normalized  by  the  total  signal  power  plus  noise  power  in  all  of  the 
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elements.  The  weight  vector,  w(y^.,r) ,  was  calculated  for  each  position  in  the  search  space,  r 

(assumed  range  r,  depth  d,  and  azimuth  0,  for  MFP  and  assumed  azimuth  0  for  PWB).  The 
correlations  were  averaged  over  the  Nf  =  10  frequencies  to  generate  the  multi-tone  correlation 
estimate 

_ _  1 

C(r,d,e,a,t)  =  —  J2Cifi,r,d,Q,a,t)  ■  (3) 

Ay  /=i 


The  adaptive  process  used  in  the  analysis  was  the  Minimum  Variance  Distortionless 
Response  (MVDR)  with  white  noise  constraint^  The  MVDR  weight  vectors  are  calculated  using 


(K+sl)'U 

w=-i - i -  , 

s*(K  +  el)-'s 


(4) 


where  K  is  the  covariance  matrix,  I  is  the  identity  matrix,  and  8  is  the  white  noise  that  is  added  to 
reduce  the  sensitivity  to  replica  mismatch"'  at  high  SNR.  The  steering  vector  s  is  given  by 


s  = 


(5) 


The  replicas,  g(/-,r) ,  are  calculated  using  a  model  of  the  acoustic  field.  For  PWB,  the  model  is 

a  plane  wave  at  frequency,  /,• ,  and  azimuthal  arrival  angle,  6.  Targets  are  located  in  the  PWB 
search  space  given  by  =  (0)  •  For  MFP,  the  replicas  are  calculated  using  Collins'  Finite 
Element  Parabolic  Equation  (FEPE)  propagation  modeP  for  the  bathymetry  of  azimuth,  0,  for 
each  frequency,  /, ,  and  for  each  range,  r,  and  depth,  d.  Targets  are  located  in  the  MFP  search 
space  given  by  =  (r,d,Q). 

The  element  positions  of  the  bottom  mounted  HLA  were  measured  prior  to  the  source 
tows.  Independent  Array  Element  Location  (AEL)  measurements  of  the  VLA  were  made  during 
the  events^  but  were  not  yet  available  for  this  analysis.  Instead  the  array  was  assumed  to  be 
straight  and  the  array  tilt,  a,  was  included  in  the  search  parameters.  The  MFP  results  presented 
here  are  for  array  tilts  corresponding  to  the  tilt  with  maximum  C{r,d, ,  t)  obtained  using  the 

§ 

pilot  tones  and  linear  (Bartlett)  MFP,  w  =  — .  Array  shapes  derived  by  the  two  methods  are 

N 

e 

compared  in  ref  6. 

The  loss  from  replica  mismatch  was  measured  using  the  pilot  tones  for  both  PWB  and 
MFP  replicas  with  the  methods  described  in  ref  2.  In  both  cases,  the  loss  was  less  than  1  dB  for 
the  ranges,  frequencies  and  arrays  used  in  the  experiment. 
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Adaptive  MFP  Results  using  VLA 


The  adaptive  MFP  beamforming  described  above  was  done  on  the  ArcMFP-5  VLA  data 
set  using  replicas  generated  with  the  bathymetry  along  the  azimuth  of  the  track  shown  in  Fig.  1. 
Displays  of  output  correlation  vs  range  and  time  were  generated  at  several  search  depths. 
Processing  was  done  using  the  ten  pilot  tone  frequencies,  the  “-22  dB”  set  and  a  noise  only  set  of 
ten  frequencies  at  fj+9  Hz,  1  Hz  above  the  “-22  dB”  set.  Fig.  2  shows  range-time  ambiguity 
surfaces  for  the  target  depth  of  30  m. 


Figure  2:  Range-Time  ambiguity  surfaces  from  VLA  for  the  pilot  tones,  the  “-22  dB”  tones  and 
the  noise  only  frequencies  at  the  target  depth  of  30  m.  MFP  using  MVDR  algorithm  with  white 
noise  constraint  with  white  noise  gain  =  6  dB.  The  VLA  consisted  of  8  elements  spaced  at 
yi  @  67  Hz. 


The  pilot  tone  surface  (far  left)  shows  the  track  of  the  target.  The  ISNR,  seen  in  Fig  1 , 
was  approximately  12  dB  for  the  first  35  minutes  of  the  track.  Sidelobes  of  the  signal  are  also 
seen.  These  sidelobes  were  not  reduced  by  the  adaptive  processing  because  of  the  white  noise 
constraint.  The  “-22  dB”  tone  surface  (center)  illustrates  the  detection  and  track  of  the  low  level 
signal  at  -  10  dB  ISNR.  The  noise  only  frequencies  surface  is  shown  at  the  right.  The  array  gain 
for  the  low  level  signal  vs  time  can  be  calculated  by 


AG  =  10  log 


/"j.) 


-ISNR 


(6) 
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C'j-(  /, 

W 

output  correlation  at  the  target  range  and  depth,  and  C^(t)is  the  average  of  the  output 

correlation  of  the  noise  only  surface  over  range  at  the  target  depth.  For  the  case  of  Fig  2,  the 
array  gain  was  the  order  of  10  dB.  This  is  consistent  with  the  “-22  dB”  track  disappearing  when 
the  ISNR  dropped  below  -10  dB. 

Examining  the  “-22  dB”  surface,  it  can  be  seen  that  the  source  track  is  discemable  in  a 
low  clutter  background.  The  low  clutter  results  from  the  depth  resolution  of  MFP  combined  with 
the  sidelobe  reduction  from  frequency  averaging  and  MVDR  processing. 


is  the  output  signal  to  noise  ratio  (OSNR),  Cy{t,rj.,dj.)  is  the 


where  10  log 


Adaptive  PWB  Results  using  HLA 

The  adaptive  PWB  beamforming  described  above  was  done  on  the  ArcMFP-5  HLA  data. 
The  HLA  was  oriented  so  that  the  broadside  direction  was  pointing  -10°  re:  true  North.  The 
PWB  replicas  provided  a  search  space  from  -90°T  to  +90°T.  Displays  of  output  correlation  vs 
azimuth  (re:  true  North)  and  time  were  generated.  Processing  was  done  using  the  ten  pilot  tone 
frequencies,  the  “-22  dB”  set  and  a  noise  only  set  of  ten  frequencies  at  yj+9  Hz,  1  Hz  above  the 
22  dB”  set.  Fig.  3  shows  azimuth-time  ambiguity  surfaces  (bearing  time  records)  for  the 
ArcMFP-5  data  set.  The  target  track  is  between  the  two  dark  lines  in  the  pilot  tone  surface. 


Figure  3:  Azimuth-Time  ambiguity  surfaces  from  HLA  for  the  pilot  tones,  the  “-22  dB”  tones 
and  the  noise  only  frequencies  at  the  target  depth  of  30  m.  PWB  using  MVDR  algorithm  with 
white  noise  constraint  with  white  noise  gain  =  6  dB.  The  HLA  consisted  of  14  elements  spaced  at 
X/2  @  187.5  Hz.  The  track  of  the  towed  source  is  between  the  dark  lines  which  were  drawn  on 
the  pilot  tone  surface. 


The  ISNR  was  measured  separately  on  the  HLA  and  was  within  1  dB  of  that  shown  in  Fig 
1  for  the  VLA,  ~  12  dB  for  the  first  35  minutes  of  the  track.  The  “-22  dB”  tone  surface  is  in  the 
center  with  an  ISNR  of-  -10  dB.  The  noise  only  frequencies  surface  is  shown  at  the  right. 

The  most  prominent  feature  of  all  surfaces  is  common-mode  electronic  noise  in  the 
broadside  direction  (-10°T).  The  common-mode  noise  was  in  the  same  direction  as  the  signal  at 
the  end  of  the  track.  This  noise  combined  with  other  effects  discussed  in  the  next  paragraph  made 
quantitative  comparison  with  the  VLA  MFP  result  using  this  data  fruitless. 

Several  interfering  ships  can  be  seen  in  the  noise  only  surface.  The  noise  only  surface  also 
shows  a  detection  track  when  the  source  is.  These  effects  illustrate  several  issues  of  importance 
in  experimentally  evaluating  the  detection  and  classification  performance  of  an  HLA  with  PWB  on 
submerged  signals. 

•  With  the  HLA,  additional  means  of  submerged  target  identification  and  classification  must  be 

used  to  distinguish  the  target  from  the  clutter. 

•  The  output  SNR  for  array  gain  estimates  should  be  taken  locally  to  exclude  interferences.  The 

noise  estimate  should  exclude  any  tow  ship  signal. 

•  The  variance  in  the  “noise  only”  search  space  relates  to  the  detection  threshold  of  the  display. 

If  a  target  is  well  separated  from  the  clutter  it  can  be  detected,  so  the  variance  should  also 

be  measured  locally. 

As  a  result  of  these  difficulties,  quantitative  comparison  is  not  definitive.  However,  it  can  be  said 
that  the  low  level  signals  were  detected  with  about  the  same  array  gain  (-10  dB)  as  the  VLA  with 
MFP. 


Issues  for  SWellEX-96 

The  SWelIEX-96  experiment  was  conducted  at  the  same  site  in  May  of  1996  to  address 
new  issues  and  to  incorporate  the  lessons  learned  in  SWellEX-3.  It  was  planned  and  conducted 
by  NRaD  in  conjunction  with  the  following  projects: 

•  Environmentally  Enhanced  Signal  Processing  (EESP)  NRL-7100 

•  All  Optical  Deployable  System  (AODS)  PMW-183 

•  USS  Dolphin  Project  PMS-395 

The  MFP  objectives  of  SWelIEX-96  are  listed  below: 

•  Demonstrate  Matched  Field  Processing  (MFP)  techniques  on  extended  targets  whose  size  is 

on  the  order  of  the  resolution  cell. 

•  Determine  the  feasibility  of  MFP  at  frequencies  up  to  700  Hz. 

•  Compare  broadband  and  multi-tone  MFP  detection  performance  between  a  tilted  vertical  line 

array  (TVLA)  and  a  vertical  line  array  (VLA).  The  TVLA  has  range,  depth  and  azimuthal 
resolution  providing  the  potential  for  adaptive  processing  to  resolve  and  null  interferences 
in  azimuth. 

•  Compare  broadband  and  multi-tone  minimum  detectable  level  of  HLA  plane  wave 

beamforming  with  VLA  and  TVLA  MFP  using  the  diesel-electric  submarine  USS 
DOLPHIN  with  augmentation. 
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The  SWellEX-96  geometry  is 
shown  in  Fig.  5.  The  acoustic 
sensors  deployed  included; 

A  vertical  line  array  (VLA)  and 
a  tilted  vertical  line  array 
(TVLA)  were  deployed  by 
FLIP,  shown  in  Fig  6. 

The  All  Optical  Deployable 
System  (AODS)  consistmg  of 
two  horizontal  line  arrays 
(HLAs)  was  installed  on  the 
sea  floor. 

The  SatelUte  linked  Vertical 
Lme  Array  (SVLA),  also 
shown  in  Fig.  6,  was  brought 
to  the  experiment  by  the  NRL 
EESP  project.  SVLA  is  a 
remotely  controlled,  self 
recording  VLA  system  which 
extends  1 65  m  up  fi'om  the  sea 
Figure  4:  SWelIEX-96  Experiment  Geometry  floor. 


Figure  5:  SWellEX-96  Vertical  Arrays  for  Matched  Field  Processing 

The  environment  was  sampled  to  support  control  and  interpretation  of  the  various 
experimental  events.  The  full  details  of  SWellEX-96  events  are  published  in  the  test  plans  and 
preliminary  data  report,  which  may  be  obtained  from  the  author.  Processing  onboard  FLIP  and  at 
the  onshore  SVLA  site  during  the  experiment  detected,  localized  and  tracked  multi-tone  signals 
with  time-late  of  3  hours,  compared  to  the  3  days  achieved  ashore  during  SWellEX-3. 
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The  difficulties  experienced  in  SWellEX-3  should  have  been  alleviated  in  the  following 

ways: 

•  The  augmented  quiet  diesel-electric  research  submarine,  USS  DOLPHIN,  carried  the  source 

transmitting  multi-tone  MDL  signals  which  spanned  frequencies  between  50  and  680  Hz. 
If  DOLPHIN'S  signature  appears  in  the  noise  only  surfaces,  it  will  appear  in  both  the  VLA 
and  HLA  surfaces  with  the  same  source  level. 

•  The  AODS  HLA  array  had  no  common-mode  electrical  noise  at  broadside  because  of  its  all 

optical  design.  The  HLAs  were  oriented  such  that  source  tracks  were  not  in  the  broadside 
direction,  so  that  any  common-mode  noise,  should  it  appear,  will  not  appear  on  the  source 
azimuths. 


Summary  and  Conclusions 

The  results  of  the  SWelIEX-3  experiment  have  demonstrated  the  potential  of  using  vertical 
line  arrays  with  adaptive  MFP  for  the  detection  and  tracking  of  submerged  signals.  Broad-band 
adaptive  processing  coupled  with  the  depth  resolution  of  MFP  shows  a  low- variance  clutter  free 
search  space  at  depth  for  signal  detection.  The  array  gain  of  the  VLA  was  comparable  to  the 
HLA  with  adaptive  PWB. 

The  SWellEX-96  experiment  was  designed  to  overcome  many  of  the  difficulties  in 
comparing  VLA  and  HLA  array  detection  performance  which  were  experienced  in  SWellEX-3. 
Several  issues  remain  to  be  addressed  to  accomplish  a  quantitative  and  objective  comparison 
including: 

•  Determining  the  probability  of  false  alarm  in  the  large  MFP  search  space. 

•  Determining  the  mean  and  variance  of  the  noise  in  the  cluttered  search  space  of  the  HLA. 

The  hypothesis  will  be  tested  that  the  range,  depth  and  azimuthal  resolution  of  the  tilted  vertical 
array  will  improve  the  detection  performance. 
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Abstract 

SWelIEx-1  (Shallow  Water  evaluation  cell  Experiment  #1)  and  SWelIEx-3  were  carried  out  in  August  1993 
and  July  1994  west  of  Point  Loma  in  approximately  200  m  water.  During  SWellEx-1,  a  MPL  48-element  vertical 
array  was  deployed  from  the  R/P  FLIP.  Similarly,  during  SWellEx-3,  a  MPL  64-element  vertical  array  was 
deployed  from  the  R/P  FLIP.  Source  tow  events  were  carried  out  throughout  the  region  and  included  both  range- 
independent  radial  tracks  as  well  as  tracks  over  a  variable  bathymetry.  In  SWellEx- 1,  a  set  of  4  tonals  were 
transmitted  in  the  50-200  Hz  band  while  in  SWellEx-3, 10  tonals  were  transmitted.  Shipping  noise  varied 
substantially  during  these  experiments  due  to  fluctuations  in  traffic  patterns. 

The  focus  of  this  paper  is  on  the  performance  of  the  MV-EPC  (Minimum  Variance  -  Environmental 
Perturbation  Constraint)  matched  field  processor  with  SWellEx-1  and  SWelIEx-3  data.  For  selected  data 
segments,  a  comparison  is  made  between  the  broadband  averaged  range-depth  ambiguity  surfaces  at  the  output  of 
the  following  processors:  (1)  conventional  (Bartlett)  matched  field  processor  given  a  measured  speed  profile  at  the 
vertical  array  and  good  estimates  of  geoacoustic  parameters,  (2)  unconstrained  MVDR  (Minimum  Variance 
Distortionless  Response)  processor  using  the  same  environmental  information  as  (1),  and  (3)  MV-EPC  processor 
given  the  measured  sound  speed  profile  at  the  vertical  array  and  a  statistical  description  of  the  geoacoustic 
environment.  Also,  range-time  ambiguity  surfaces  (at  the  depth  of  the  source)  are  shown  to  illustrate  the  time- 
evolving  structure  of  the  surfaces  at  the  output  of  these  processors. 

pVork  supported  by  ONR,  Code  321US.] 
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Processors 


•  Bartlett 

•  MVDR  (Minimum  Variance  Distortionless  Response) 

•  MVDR  -  NLC  (Neighborhood  Location  Constraints) 

•  MVDR  -  WNC  (White  Noise  Constraint) 

•  MVDR  -  EPC  (Environmental  Perturbation  Constraints) 
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The  SWeIlEx-1  environment. 


Environmental  Parameter 

Lower 

Bound 

Baseline 

Value 

Upper 

Bound 

Units 

Bottom  depth 

195.0 

197.6 

200.0 

m 

Sediment  velocity  (top) 

1545.0 

1554 

1565.0 

m/s 

Sediment  velocity  (bottom) 

1565.0 

1576 

1585.0 

m/s 

Sediment  density 

1.65 

1.74 

1.85 

g/cm^ 

Sediment  attenuation 

0.1 

0.2 

0.3 

dB/(km  Hz) 

Sediment  thickness 

28.5 

30.0 

32.5 

m 

Bottom  velocity 

1850.0 

1861.0 

! 

1870.0 

m/s 

Bottom  density 

1.9 

2.0 

2.1 

Bottom  attenuation 

0.1 

1 

0.20 

0.3 

dB/(km  Hz) 

SWelIEx-1  Environmental  Parameters. 
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SWellEx-1  -  Source  Parameters 


•  Depth  =  82  m 

•  Speed  =  3  knots 

•  Frequencies  =  70,  95,  145,  and  195  Hz 


83 


Lofargram  of  the  G-I  Track  event 
(source  frequencies  of  70, 95, 145,  and  195  Hz). 


Z^bart^favg.mat.  SD  =  82  m  —  JD  229,  1 6: 1 5-1 7:36  PDT  (Max  =  0  dB) 


^ISIHBRasi 


Range  (m) 

Z_mv_favg.mat.  SD  =  82  m  —  JD  229,  16:15-17:36  PDT  (Max  ^ 


8 

Range  (m) 

Z_mvepc_favg.mat,  SD  =  82  m  —  JD  229,  16:15-17:36  PDT  (Max 

=  -7 

dB) 

Magnitude  (dB) 

Time-range  ambiguity  surfaces  for  the  Bartlett,  MV,  and  MV-EPC  methods 

(frequency-averaged). 
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SWellEX-3  ArcMFP-5  Event,  94.Jul.30 

DGPS  track  of  R/V  Acoustic  Explorer,  211.1300Z-211.1540Z 


isobath  intervals  20  znetera 


chart  pr^auccd  hj  KCCOSC  KOTIcE  Mt.  Code  641 


SWellEx-3  Vertical  Array 


Transponder 


Depth  (Meiers) 


SWellEX-3  CTD  Sound  Speed  Profile  No.  151,  Site  G,  30  July  1994,  1319  ZULU 


1485  1490  1495  1500  1505  1510  1515  1520  1525 


1485  1490  1495  1500  1505  1510  1515  1520  1525 

Sound  Speed  (M/Sec) 
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The  SWeIlEx-3  environment. 
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Environmental  Parameter 

Lower 

Bound 

Baseline 

Value 

Upper 

Bound 

Units 

Bottom  depth 

195.0 

198.0 

200.0 

m 

Sediment  velocity  (top) 

1560.0 

1572.368 

1580.0 

m/s 

Sediment  velocity  (bottom) 

1585.0 

1593.016 

1605.0 

m/s 

Sediment  density 

1.65 

1.76 

1.85 

g/cm^ 

Sediment  attenuation 

0.1 

0.2 

0.3 

dB/(kmHz) 

Sediment  thickness 

28.5 

30.0 

32.5 

m 

Mudstone  velocity  (top) 

1870.0 

1881.0 

1890.0 

m/s 

Mudstone  velocity  (bottom) 

3235.0 

3245.8 

3255.0 

m/s 

Mudstone  density 

1.95 

2.06 

2.15 

g/cm^ 

Mudstone  attenuation 

0.01 

0.06 

0.11 

dB/(km  Hz) 

Mudstone  thickness 

798.0 

800.0 

802.0 

m 

Basement  velocity 

5200.0 

5200.0 

5200.0 

m/s 

Basement  density 

2.66 

2.66 

2.66 

g/cm^ 

Basement  attenuation 

0.02 

0.02 

0.02 

dB/(km  Hz) 

SWellEx-3  Environmental  Parameters. 
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SWellEx-3  -  Source  Parameters 


•  Depth  =  30  m 

•  Speed  =  5  knots 

•  Frequencies  =  53,  69,  85,  101,  117,  133,  149,  165,  181,  and 
197  Hz 


93 


Magnitude  (dB) 


ESR:  3 100  m,  ESD;  30  m  (Max  =  -27  clB)  Z_mvnlc  Javg.8  ~  ESR 


Magnitude  (dB) 


,wnc_favg.O  —  ESR:  3060  m,  ESD:  30  m  (Max  =  -39  dB)  Z_wnc_favg.8  —  ESR:  4140  m,  ESD:  30  m  (Max  =  -39 
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Matched  Field  Inversion  of  Real  Acoustic 
Data  Using  Simulated  Annealing  and 
Genetic  Algorithms 
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Abstract 

Localization  of  underwater  acoustic  sources  by  matched  field  processing  requires  very  accurate 
knowledge  of  the  ocean  environment.  This  stringent  requirement  is  circumvented  in  the  localization 
approach,  in  which  environmental  parameters  are  included  in  the  parameter  search  space.  Due 
to  the  extremely  large  number  of  parameter  combinations,  this  inverse  problem  requires  modem 
global  optimization  methods. 

Experimental  results  on  the  simultaneous  estimation  of  source  position,  geo-acoustic  parameters 
and  other  parameters  (such  as  water  depth  and  array  geometry)  are  presented  for  a  shallow  water 
environment.  The  stationary  source  at  a  range  of  about  5.4  km  from  the  vertical  array  is  succesfully 
localized  in  range  and  depth. 

For  the  forward  propagation  model  we  used  a  normal-mode  code  developed  at  TNO-FEL.  For 
the  global  optimization  both  simulated  annealing  algorithms  and  genetic  algorithms  have  been 
employed.  The  performance  of  these  algorithms  with  regard  to  efficiency  and  robustness  have 
been  investigated.  Further,  different  parametrizations  of  the  ocean  environment  (with  different 
degrees  of  complexity)  have  been  employed. 

The  inversion  is  followed  by  a  statistical  analysis  of  the  parameter  values  obtained.  When  simulated 
annealing  is  used  for  the  optimization  this  has  been  accomplished  by  running  the  algorithm  a 
large  number  of  times.  Both  the  uncertainty  in  the  inverted  parameters  (which  is  related  to  their 
acoustical  importance)  and  the  correlation  between  the  inverted  parameters  have  been  investigated. 

1  Introduction 

During  the  last  decade  Matched  field  processing  (MFP)  has  become  an  important 
research  topic  in  underwater  acoustics.  The  main  application  of  MFP  is  localization 
of  underwater  acoustic  sources. 

The  basic  principle  of  MFP  is  simple:  The  acoustic  field  measured  with  a  (vertical) 
hydrophone  array  is  correlated  with  the  acoustic  field  computed  by  an  appropri¬ 
ate  propagation  model  for  candidate  values  of  the  range  and  depth  of  the  source. 
This  yields  the  so-called  range/depth  ambiguity  surface.  The  candidate  range  and 
depth  showing  the  highest  correlation  between  measured  and  modelled  field  should 
correspond  to  the  true  range  and  depth  of  the  source. 

In  MFP  it  is  assumed  that  all  environmental  (and  other)  input  parameters  for  the 
model  are  accurately  known.  However,  these  parameters  are  generally  not  known  or 
at  least  not  sufficiently  accurate  for  a  successful  source  localization.  This  situation 
is  referred  to  as  mismatch. 
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Under  practical  conditions,  i.e.  when  using  real  experimental  acoustic  data  obtained 
at  sea,  one  always  encounters  mismatch.  A  very  promising  method  for  solving  this 
problem  and  for  making  MFP  more  robust  (a  requirement  for  practical  applications) 
is  the  so-called  focalization  approach  [1].  In  this  case  both  the  source  position  and 
the  environmental  parameters  are  candidate  parameters,  i.e.  a  range,  depth  and 
environment  are  sought  for  which  the  acoustic  field  focuses  into  a  source  point. 
Then,  environmental  parameters  and  other  (geometrical)  parameters  such  as  array 
position,  are  being  estimated  as  part  of  the  localization  process.  Due  to  the  ex¬ 
tremely  large  number  of  possible  parameter  combinations,  such  an  inverse  problem 
requires  modern  global  optimization  methods  such  as  simulated  annealing  [1]  or 
genetic  algorithms  [4]. 

Hereafter  focalization  is  referred  to  as  matched  field  inversion. 

The  success  of  MFP  stems  also  from  the  availability  of  highly  accurate  propaga¬ 
tion  models.  At  present,  normal-mode  models  are  the  most  appropriate  for  MFP 
calculations  [7]. 

In  this  paper  results  are  presented  for  real  experimental  data  obtained  in  shallow 
water.  For  the  propagation  model  we  have  used  a  normal-mode  model  developed  at 
TNO-FEL  [6],  whereas  both  simulated  annealing  and  genetic  algorithms  have  been 
used  as  a  global  optimization  method. 

Section  2  provides  a  brief  description  of  the  experimental  setup,  i.e.  the  environ¬ 
mental  conditions  of  the  trial  area  and  the  characteristics  of  the  transmitting  and 
receiving  equipment. 

Section  3  discusses  the  basic  principle  of  simulated  annealing  (SA)  and  the  results 
obtained  for  the  parameter  inversion,  whereas  section  4  discusses  the  basic  principle 
of  genetic  algorithms  (GA)  and  the  results  obtained  with  this  technique. 

In  section  5  the  conclusions  of  this  work  are  given. 

2  Experimental  setup 

The  experimental  data  we  used  were  obtained  from  SACLANTCEN.  In  October 
1993  SACLANTCEN  has  performed  an  MFP  experiment  in  a  shallow  water  area 
north  of  the  island  of  Elba  (position  43°02.86'N,  10°10.01'E).  In  this  area  the  geo¬ 
acoustic  properties  were  reasonably  well  known  from  previous  experiments,  i.e.  a 
priori  knowledge  of  the  acoustic  bottom  parameters  (hereafter  denoted  by  historical 
parameters)  was  available.  The  experiment  is  described  in  great  detail  in  [5]  and  is 
briefly  summarized  here. 

The  experiment  has  been  carried  out  in  a  virtually  range-independent  area.  The 
water  depth  amounts  to  127  m  and  the  flat  bottom  is  covered  with  a  thin 
clay/sand-clay  sediment  of  a  few  meters  thickness.  The  propagation  conditions  are 
typical  for  the  summer  conditions  encountered  in  the  area  around  Elba  (isothermal 
upper  layer  of  ~  60  m).  For  the  part  of  the  data  we  selected,  the  wind  speed 
amounts  to  ~  8  m/s  corresponding  to  a  rms  surface  roughness  of  about  0.5  m. 

The  complete  sound  speed  profile  of  the  ocean  environment  is  given  in  figure  1. 
The  sound  speed  profile  in  the  water  column  has  been  measured  at  the  site  of  the 
receiving  array. 

Here  we  have  assumed  that  the  bottom  below  the  sediment  (denoted  as  subbot¬ 
tom)  is  a  semi-infinite  homogeneous  medium.  This  is  in  accordance  with  the  usual 
normal-mode  modelling  of  the  sea  floor:  the  ocean  bottom  is  assumed  to  consist  of  a 
sediment  layer  over  a  semi-infinite  homogeneous  subbottom.  In  the  sediment  layer 
the  sound  speed  is  allowed  to  vary  with  depth,  whereas  in  the  subbottom  the  sound 
speed  Cb  is  assumed  to  be  constant.  The  density  ps  and  attenuation  coefficient  as 
of  the  sediment  are  also  assumed  to  be  constant.  The  same  holds  for  pb  and  a^,  the 
subbottom  density  and  attenuation  coefficient,  respectively. 
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Figure  1:  The  ocean  environment  discussed  in  this  paper. 


Previous  measurements  carried  out  in  this  area  have  indicated  that  the  sound  speed 
in  the  sediment  increases  linearly  with  depth.  Hence,  the  sound  speed  profile  in  the 
sediment  can  be  characterized  by  three  parameters:  the  upper  sediment  speed  Cgu, 
the  lower  sediment  speed  Cgi  and  the  sediment  thickness  Hs-  Historical  values  for 
the  geo-acoustic  parameters  are  taken  from  [5]  and  are  given  in  table  I  below.  These 
values  were  also  used  in  figure  1. 

TABLE  I:  Historical  geo- acoustic  parameters  and  nominal  geometrical  parameters 


Parameter 

symbol 

rmit 

numerical  value 

upper  sediment  speed 

^su 

m/s 

1520 

lower  sediment  speed 

Ciu 

m/s 

1580 

sediment  thickness 

Hs 

m 

2.5 

sediment  density 

Ps 

g/cc 

1.7 

subbottom  density 

Pb 

g/cc 

1.8 

sediment  attenuation 

dig 

dB/A 

subbottom  attenuation 

O^b 

dB/A 

subbottom  speed 

Cb 

m/s 

source  range 

rs 

km 

5.4 

source  depth 

Zs 

m 

75 

water  depth 

Hy, 

m 

127 

array  depth 

dR 

m 

112.7 

The  receiving  system  is  a  vertical  array  consisting  of  48  hydrophones  with  2  m  spac¬ 
ing.  The  nominal  depth  of  the  deepest  hydrophone,  hereafter  denoted  as  receiving 
array  depth  dn,  is  112.7  m.  Due  to  array  tilt  and/or  imprecise  measurement  of  the 
water  depth  the  actual  array  depth  can  be  different. 
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An  acoustic  source  was  deployed  at  a  range  of  5.4  km  north  of  the  receiving 
array  at  a  depth  25  of  ~  75  m.  Two  broadband  signals  were  transmitted  by  the 
source,  not  simultaneously.  The  frequency  bands  of  these  signals  were  160  -  180 
Hz  and  320  -  350  Hz,  respectively.  For  the  analysis  reported  hereafter  we  have 
selected  a  single  frequency  in  the  centre  of  the  160  -  180  Hz  band,  i.e.  the  Fourier 
coefScients  at  168.95  Hz  obtained  after  a  Fourier  transformation  of  8  seconds  of  the 
data  received  at  each  hydrophone  (sample  frequency  of  the  data  acquisition  system 
is  1  kHz).  This  measured  complex  pressure  vector  has  been  inverted  using  simulated 
annealing  (section  3)  and  genetic  algorithms  (section  4), 

3  Inversion  by  simulated  annealing 

Matched  field  inversion  has  been  applied  to  the  measured  complex  pressure  field 
described  in  the  former  section.  Both  source  range  and  depth  and  other  geometrical 
and  environmental  model  input  parameters  are  considered  unknown  parameters, 
which  have  to  be  determined  through  an  optimization  procedure.  This  involves 
finding  a  set  of  parameter  values  which  minimizes  the  discrepancy  between  the 
modelled  and  measured  acoustic  field. 

As  a  start  the  optimization  was  carried  out  for  all  12  model  input  parameters  given 
in  table  1.  There  was  no  optimization  for  the  sound  speed  profile  in  the  water 
column:  the  profile  measured  at  the  array  site  was  used. 

The  determination  of  this  relatively  large  number  of  unknown  input  parameters  is 
a  real  inverse  problem.  The  number  of  parameter  value  combinations  is  extremely 
large.  In  addition,  the  parameter  search  space  can  have  a  large  number  of  local 
minima.  Such  a  problem  requires  modern  global  optimization  methods. 

The  application  of  simulated  annealing  as  a  global  optimization  method  for  solving 
the  inverse  problem  is  discussed  in  this  section,  whereas  the  application  of  a  genetic 
algorithm  is  dicussed  in  the  section  thereafter.  Both  techniques  are  efficient  methods 
for  finding  the  global  minimum  of  a  cost  or  energy  function  that  depends  on  many 
parameters. 

3.1  Basic  principle 

In  the  SA  algorithm  the  unknown  parameters  have  been  determined  by  minimizing 
the  energy  function  (or  cost  function) 

(1) 

with  Piin  the  linear  processor  power  given  by 

Plin  =  Weal  (m)  ■  (2) 

i.e.  the  inner  product  of  the  calculated  pressure  field  Pcai  suid  the  observed  pressure 
field  Po6s5  both  in  the  frequency  domain.  Both  column  vectors  are  normalized  to 
unity  and  are  of  length  48  (being  the  number  of  hydrophones) 


Peal  —  1  (3) 

Peal  has  been  computed  by  a  normal-mode  model  [6]  and  depends  on  m,  the  input 
parameter  vector  given  by 

771  =  (c^n ,  Cgi ,  Hs 'i  Ps')  Pb-i  5  ^6 Pw’)  )  (4) 

The  search  is  initiated  with  random  values  for  the  parameters  rrii  within  the  imposed 
bounds  The  search  bounds  for  each  model  input  parameter,  i.e.  the 

priori  knowledge,  are  given  in  table  II  below. 
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TABLE  II:  Model  input  parameters,  historical  values,  search  hounds  and  SA  results 


wmi 

By 

SA  (log) 

SA  (lin) 

Csu  (m/s) 

■EMjl 

1550 

1487  ±25 

1492±24 

Csi  (m/s) 

1600 

1535  ±25 

1541±21 

Hs  (m) 

2.5 

1.5 

6.0 

3.9  ±1.2 

4.0±1.3 

Ps  (g/cc) 

1.7 

1.2 

2.2 

1.9  ±0.3 

1.9±0.2 

Pb  (g/cc) 

1.8 

1.2 

2.2 

1.6  ±0.3 

1.8±0.3 

a,  (dB/A) 

0.0 

0.4 

0.04  ±0.05 

0.04±0.04 

at,  (dB/A) 

0.0 

0.4 

0.13  ±0.11 

0.19±0.11 

Cb  (m/s) 

1550 

1650 

1577  ±5 

1577±6 

rs  (m) 

1000 

11000 

5362  ±99 

5361±65 

Zs  (m) 

75 

10 

no 

75.5  ±0.7 

75.5±0.8 

Hy,  (m) 

125 

130 

127.5±1.6 

127.3±1.3 

dR  (m) 

112.7 

110 

114 

111.9±0.7 

111.7±0.5 

The  energy  function  is  evaluated  for  the  initial  parameter  combination. 

Then  the  individual  parameters  are  randomly  perturbed  one  at  a  time  according  to 

m[  +  ^  Ai  (5) 

with  ^  a  random  number  drawn  from  a  uniform  distribution  on  [“!,!]. 

The  maximum  perturbation  Aj  allowed  for  parameter  has  been  taken  half  the 
corresponding  search  interval 

Ai  =  -  Bi,i)  (6) 

When  the  new  parameter  value  m[  falls  outside  the  interval  a  new  value 

for  m[  is  drawn  again. 

After  each  perturbation  the  energy  function  is  evaluated.  A  decrease  in  E  {AE  <  0) 
is  accepted  unconditionally,  whereas  an  increase  in  E  {AE  >  0)  is  accepted  with  a 
probability  P  given  by  the  Boltzmann  distribution 

P  =  (7) 

with  Tj  a  control  parameter  analogous  to  the  temperature. 

When  all  parameters  have  been  perturbed  a  certain  number 
this  work  [2]),  the  temperature  is  slightly  reduced  according  to 
Two  different  sequences  of  decreasing  temperatures  have  been 
logarithmic  schedule 

Tj-f  1  =  Cf  Tj  (8) 

and  the  linear  schedule 

Tj  =  Ti  +  a  0*  -  1)  (Tend  -  Ti)  (9) 

The  cooling  factor  c/  for  logarithmic  cooling  is  a  constant  <  1,  but  close  to  1. 
For  linear  cooling  Ti  and  Tend  are  the  starting  temperature  and  end  temperature, 
respectively,  and  a  is  another  constant. 

Hence,  the  probability  of  accepting  an  increase  in  E  decreases  as  T  is  lowered  (see 
equation  7).  The  concept  of  accepting  perturbations  that  increase  E  allows  the 
algorithm  to  escape  from  local  minima  in  the  parameter  search  space. 


of  times  (twice  in 
a  cooling  schedule, 
employed,  i.e.  the 
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Figure  2:  Convergence  of  one  of  the  SA  runs. 


SA  is  problem-specific:  appropriate  values  for  Ti  and  c/  (in  the  case  of  logarithmic 
cooling)  and  for  Ti ,  Tend  a-nd  a  (for  linear  cooling)  have  to  be  determined  by  trial  and 
error.  The  cooling  must  be  sufficiently  slow  in  order  to  reach  the  global  minimum, 
thereby  avoiding  sub-optimal  solutions  (local  minima) . 

3.2  Results 

The  SA  algorithm  has  been  run  25  times  for  both  cooling  laws  mentioned  above. 
For  each  run  the  same  8  seconds  of  experimental  data  was  used. 

For  logarithmic  cooling  the  SA  parameters  c/  and  Ti  were  0.97  and  0.06,  respec¬ 
tively.  For  linear  cooling  the  SA  parameters  Ti,  Tend  and  a  were  0.06,  0.0003 
and  0.00503,  respectively.  With  400  iterations  (see  below)  then  the  start  and  end 
temperatures  of  both  cooling  sequences  are  equal. 

The  starting  temperature  has  been  determined  such  that  about  85  %  of  the  per¬ 
turbations  are  accepted  at  the  first  iteration,  see  [2].  About  200  temperature  steps 
turned  out  to  be  sufficient  for  converging  of  the  algorithm.  Since  each  parameter 
is  changed  twice  at  each  temperature,  this  corresponds  to  400  iterations.  (After  all 
parameters  have  been  changed  once,  an  iteration  is  completed). 

The  convergence  of  the  SA  algorithm  is  illustrated  in  figure  2.  In  this  figure  the  mean 
energy  averaged  over  all  accepted  model  runs  at  a  single  temperature  is  plotted  as  a 
function  of  iteration  for  one  of  the  SA  runs  (logarithmic  cooling).  Note  that  in  this 
figure  1  —  Piin  has  been  plotted  instead  of  the  (mean)  energy  1  —  (see  equation 
1)- 

The  convergence  of  the  model  input  parameters  is  given  in  figure  3.  One  observes 
that  both  the  source  range  and  source  depth  converge  rapidly  to  correct  values. 
Employing  the  logarithmic  cooling  schedule  (with  Cf  =  0.97)  the  SA  algorithm  has 
been  successful  15  times,  i.e.  10  times  the  algorithm  converged  to  a  local  minimum. 
In  the  latter  10  cases  correct  values  for  source  range  and  depth  were  not  found. 
Employing  the  linear  cooling  schedule  the  number  of  successful  runs  amounts  to 
23  (out  of  25).  Since  the  start  and  end  temperature  of  the  cooling  schedules  are 
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Figure  3:  Convergence  of  the  model  input  parameters  for  the  SA  run  of  figure  2. 


equal,  the  much  better  performance  of  the  linear  cooling  schedule  must  be  due  to 
the  slower  cooling  in  the  first  part  of  the  linear  cooling  schedule. 

In  this  way  15  estimates  for  each  model  input  parameter  are  produced  in  the  case  of 
the  logarithmic  cooling  schedule  and  23  estimates  for  each  parameter  are  produced 
in  the  case  of  the  linear  cooling  schedule.  The  a  posteriori  probability  distributions 
of  the  parameter  estimates  are  given  in  figure  4  (using  the  linear  cooling  schedule 
results).  In  this  way  information  is  provided  about  the  uncertainty  in  the  final  re¬ 
sults  obtained.  Table  II  provides  the  mean  and  standard  deviation  of  the  parameter 
estimates. 

It  is  seen  that  the  source  range  the  source  depth  and  the  subbottom  sound 
speed  C6  are  quite  well  determined.  For  these  parameters  the  a  posteriori  probability 
distributions  are  unambigously  peaked  and  the  standard  deviation  (see  table  II)  is 
small  compared  to  the  corresponding  search  intervals  (the  a  priori  information).  The 
parameters  '2^5  and  ct  turn  out  to  be  the  acoustically  most  important  parameters. 
The  other  parameters  are  less  well  determined.  Since  the  sediment  thickness  is 
less  than  half  a  wavelength  at  the  selected  frequency,  propagation  is  less  sensitive 
to  the  sediment  parameters.  The  corresponding  probability  distributions  are  not 
unambigously  peaked  (see  figure  4)  and  the  standard  deviations  turn  out  to  be 
a  substantial  fraction  of  the  corresponding  search  intervals  (see  table  II).  Conse¬ 
quently,  the  sediment  parameters  are  not  well  determined  at  the  test  frequency. 
The  same  holds  for  the  subbottom  density  and  attenuation  constant.  However,  the 
subbottom  sound  speed  Cb  turns  out  to  be  a  very  important  acoustic  parameter. 
This  is  not  surprising,  since  the  number  of  normal  modes  is  highly  dependent  on  Cb- 
It  is  also  worthwhile  examining  the  correlation  between  the  model  input  parameter 
estimates.  It  is  observed  that  the  geometrical  parameters  (r^,  2^5,  Hyj  and  dR)  are 
highly  correlated.  For  instance,  the  highest  correlation  is  observed  between  source 
range  Vg  and  water  depth  A  plot  of  the  source  range  estimates  against  water 
depth  estimates,  see  figure  5,  indeed  reveals  the  high  correlation  between  these 
parameters.  This  means  that  there  exists  a  long  valley  in  the  parameter  search 
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Figure  5:  Source  range  estimates  ploted  against  water  depth  estimates  obtained  from 
the  SA  runs  with  linear  cooling. 


space  that  is  oriented  obliquely  to  the  parameter  axes.  Such  a  feature  is  difficult  to 
navigate.  However,  rotating  the  the  coordinate  system  such  that  the  valley  aligns 
with  one  of  the  axes  can  improve  the  efficiency  of  the  SA  algorithm  [3]. 

As  mentioned  above,  propagation  cannot  be  very  sensitive  to  the  sediment  param¬ 
eters  at  the  present  test  frequency.  Therefore,  the  SA  algorithm  has  also  been  run 
using  only  7  input  parameters.  In  this  case  the  sediment  parameters  are  omitted, 
thereby  decreasing  the  complexity.  At  the  same  time  the  parameter  search  intervals 
were  enlarged. 

Using  the  logarithmic  cooling  law  (with  c/  =  0.97)  the  algorithms  has  been  sucessful 
for  25  times  out  of  45  runs  (success  rate  62  %).  In  these  cases  both  source  range  and 
source  depth  converge  rapidly  to  correct  values,  whereas  the  third  most  important 
parameter,  subbottom  sound  speed,  converges  to  a  value  in  agreement  with  that 
obtained  for  the  SA  runs  with  12  parameters. 

The  7  model  input  parameters,  the  corresponding  search  bounds  imposed  and  the 
means  and  standard  deviations  of  the  estimated  parameter  values  obtained  after 
the  25  successful  SA  runs  are  given  in  table  III  below. 

TABLE  III:  Model  input  parameters,  search  bounds  and  estimates  resulting  from 
the  SA  runs  with  only  7  parameters 


TUi 

Bt 

Bu 

SA  estimates 

Pb  (g/cc) 

2.6 

2.3±0.3 

a6  (dB/A) 

1.0 

0.21±0.11 

Cb  (m/s) 

1750 

1575±5 

rs  (m) 

0 

11000 

5435±75 

Zs  (m) 

0 

120 

75.9±0.7 

(m) 

125 

135 

129.6±1.3 

dR  (m) 

110 

114 

112.2±0.9 
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Figure  6:  Probability  distributions  of  the  parameter  estimates  of  the  SA  runs  with 
only  7  parameters. 


The  a  posteriori  probability  distributions  of  the  parameter  estimates  are  given  in 
figure  6.  These  results  are  in  accordance  with  the  results  obtained  from  the  SA  runs 
with  12  parameters  (see  table  II). 

It  is  interesting  to  note  that  now  all  parameters  are  highly  correlated  with  each 
other  (except  for  the  bottom  attenuation  constant). 

4  Inversion  by  a  Genetic  Algorithm 

4.1  Basic  principle 

The  genetic  Algorithm  (GA)  [4]  starts  with  a,  randomly  generated,  initial  pop¬ 
ulation  of  q  members,  each  representing  a  possible  parameter  value  combination 
m.  Depending  on  the  usefulness  of  the  solution  represented,  all  members  of  the 
population  are  assigned  a  fitness  score.  Based  on  these  fitness  scores,  a  parental 
distribution  is  selected  from  the  initial  population.  The  parents  are  combined  in 
pairs,  and  operators  are  applied  to  them  to  produce  children,  i.e.  a  new  set  of  pos¬ 
sible  solutions.  Traditionally,  the  operators  are  crossover  and  mutation.  In  order  to 
apply  these  operators  the  parents  are  coded  in  binary  format.  The  parents  in  binary 
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form  axe  called  chromosomes,  consisting  of  genes,  with  each  gene  representing  one 
parameter.  If  now  a  fraction  of  the  current  population  is  replaced  by  the  children, 
a  new  population  is  created  that  contains  a  higher  proportion  of  the  characteristics 
possessed  by  the  good  members  of  the  previous  generation. 

By  continuing  this  process  over  many  generations,  good  characteristics  are  spread 
throughout  the  population,  meanwhile  being  mixed  and  exchanged  with  other  good 
characteristics. 

For  each  of  the  members  of  the  population  the  cost  or  energy  function 

E  =  (10) 

is  calculated  (the  fitness  equals  (1  —  E)),  Pun  is  given  by  equation  2.  The  parental 
distribution  is  selected  by  assigning  a  selection  probability  to  each  of  the  members 
according  to 


exp(-^^) 

T,l=iexpi-^) 


(11) 


The  temperature  T  is  decreased  as  the  evolution  continues. 

In  the  reproduction  part  q  children  are  produced.  The  first  part  of  the  reproduc¬ 
tion  is  crossover.  For  each  pair  of  parents  crossover  points  at  each  of  the  genes 
in  a  chromosome  are  randomly  selected,  and  the  gene  parts  of  both  parents  are 
exchanged.  Crossover  occurs  with  a  probability  pc-  Following  crossover  each  bit  of 
the  chromosome  is  perturbed  with  a  mutation  probability  Pm- 
Out  of  the  set  of  q  children,  (1  —  f)q  children  are  randomly  selected  to  replace  the 
{1  —  f)q  least  fit  members  of  the  initial  population.  Hence,  the  fq  most  fit  members 
of  the  initial  population  are  retained. 

This  process  of  producing  a  new  generation  is  continued  until  a  specified  amount 
of  generations  is  reached. 


4.2  Results 

The  GA  has  been  applied  to  the  inversion  process  using  the  same  data  as  used  in 
SA.  Following  [4]  the  GA  has  been  set  as  follows: 

•  The  population  size  q  should  be  large  enough  that  the  model  parameter  vectors 
can  represent  sufficient  points  in  the  parameter  search  space,  but  it  should  be 
small  enough  that  several  iterations  can  be  performed.  A  population  size  of 
64  members  is  used; 

•  The  reproduction  size  /  should  be  large  enough  such  that  the  fittest  members 
remain  in  the  population,  but  it  should  also  be  small  enough  in  order  to  avoid 
local  minima.  Here  two  different  values  for  /  have  been  examined,  i.e.  /  =  0,5 
[4]  and  /  =  0.2; 

•  In  each  generation  the  temperature  is  taken  equal  to  the  lowest  of  the  energy 
function  values  found  (i.e.  highest  fitness).  Hence,  the  temperature  decreases 
as  the  evolution  continues,  which  results  in  a  stricter  reproduction  selection 
criterion  (i.e.  a  pi  vs.  E  distribution  more  peaked  towards  higher  E,  see 
equation  11); 

•  For  the  crossover  probability  Pc  a  value  of  0.8  is  used; 

•  A  value  of  0.05  is  taken  for  the  mutation  probability 

•  The  number  of  generations  used  is  80,  which  is  sufficient  for  converging  of  the 
algorithm; 
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Figure  7:  Convergence  of  one  of  the  GA  inversions. 


•  All  parameters  are  represented  by  an  8  bit  binary  string. 

The  genetic  algorithm  has  been  run  15  times  for  /  =  0.5,  Only  40  %  of  these 
runs  converged  to  the  correct  values  for  source  range  and  the  depth.  For  /  = 
0.2,  however,  75%  of  the  runs  were  successful.  Therefore,  this  situation  has  been 
investigated  more  elaborately  (55  runs).  The  better  performance  of  the  algorithm 
with  the  lower  /  is  due  to  the  fact  that  it  allows  a  more  extensive  exploration  of 
the  parameter  search  space  (thereby  avoiding  local  minima). 

The  convergence  of  the  GA  algorithm  is  illustrated  in  figure  7.  In  this  figure  the  en- 
ergy  of  the  most  fit  member  of  the  population  is  plotted  as  a  function  of  generation. 
Note  that  1  —  Pun  is  plotted  instead  of  energy  \/l  -  Pun^ 

The  corresponding  convergence  of  the  model  input  parameters  of  this  GA  run  is 
given  in  figure  8. 

In  the  final  population  about  30  %  of  the  members  have  converged  (i.e.  (1  —  Pun)  < 
0.07).  The  probability  distributions  of  these  most  fit  members  are  given  in  figure 
9.  This  result  cannot  be  compared  directly  to  that  obtained  with  the  SA  runs, 
see  figure  6,  since  the  results  of  several  GA  runs  are  required  in  order  to  obtain 
convergence  to  other  optimum  parameter  combinations.  Still,  a  single  GA  run 
provides  information  about  the  uncertainty  (acoustical  importance)  of  the  inverted 
parameters. 

5  Conclusions 

It  has  been  demonstrated  that  both  simulated  annealing  and  genetic  algorithms  are 
capable  of  inverting  shallow  water  acoustic  data.  In  particular  the  robustness  of  the 
algorithms,  i.e.  the  relative  number  of  successful  inversions,  has  been  investigated. 
The  parameters  of  the  algorithms  have  been  established  such  that  the  success  rate 
is,  more  or  less,  maximized  without  losing  too  much  speed  of  convergence. 
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Figure  8:  Convergence  of  the  model  input  parameters  for  the  GA  run  shown  in 
figure  7. 


It  has  been  emphasized  that  an  inversion  is  incomplete  if  not  followed  by  a  statistical 
analysis  of  the  model  input  parameter  estimates.  Such  an  a  posteriori  analysis  pro¬ 
vides  information  about  the  uniqueness  and  uncertainty  (i.e.  acoustical  importance) 
of  the  parameter  estimates.  It  is  also  useful  to  investigate  parameter  coupling,  i.e. 
correlation  between  the  inverted  parameters. 
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Figure  9:  Probability  distributions  of  the  most  fit  members  of  the  final  population 
(GA  run  of  figure  7) 
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Abstract 

In  many  estimation  problems,  the  set  of  unknown  parameters  can  be  divided  into  a  subset  of  desired 
parameters  and  a  subset  of  nuisance  parameters.  Using  a  maximum  a  posteriori  (MAP)  approach  to 
parameter  estimation,  these  nuisance  parameters  are  integrated  out  in  the  estimation  process.  This 
can  result  in  an  extremely  computation  ally-intensive  estimator.  In  this  paper,  we  present  the  fast, 
approximate  MAP  parameter  estimator,  FASTMAP,  As  an  example,  we  apply  the  fast  algorithm  to 
matched-field  source  localization  in  an  uncertain  environment- 


Introduction 

In  many  estimation  problems,  the  set  of  unknown  parameters  can  be  divided  into  two  subsets.  The  param¬ 
eters  of  major  interest  comprise  one  subset  and  the  remaining  parameters  comprise  the  second  subset.  The 
remaining  parameters  only  serve  to  complicate  the  problem  and  are  referred  to  as  nuisance  parameters.  In 
the  maximum  o  posteriori  (MAP)  approach  to  parameter  estimation,  these  nuisance  parameters  are  treated 
as  random  variables  with  assumed  prior  probability  density  functions  and  integrated  out  in  the  process  of 
estimating  the  desired  parameters  [1],  [2].  Depending  on  the  number  of  nuisance  parameters,  practical  im¬ 
plementation  of  the  MAP  estimator  can  be  extremely  computationally-intensive.  In  the  next  section,  we  will 
show  that  under  certain  conditions  a  computationally-efficient  approximation  to  the  MAP  estimator  can  be 
obtained.  Efficiency  is  achieved  by  approximately  performing  the  integration  off-line  prior  to  the  processing 
of  data  observations. 

Now  let  us  describe  the  MAP  estimator  following  the  development  of  Richardson  and  Nolte  [3].  Consider 
the  following  data  model 

y  =  sa(0,^)+n,  (1) 

where  y  is  a  A'  x  1  vector  of  observed  data,  s  is  a  Gaussian  distributed  complex  random  variable,  a(0,  ^) 
is  a  X  1  vector  parameterized  by  the  vectors  0  and  and  n  is  a  A'  x  1  vector  of  Gaussian  random 
variables.  This  model  is  applicable  to  array  processing  and  time-series  analysis  problems.  We  will  assume 
that  the  parameters  of  interest  are  contained  in  0,  while  ^  contains  the  nuisance  parameters.  Eaich  of  the 
nuisance  parameters  is  assumed  to  be  a  random  variable  with  a  known  uniform  probability  distribution.  We 
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can  then  express  the  a  posteriori  probability  density  function  (pdf)  of  0  as 


p(®ly)  =  (2) 

The  data  vector  y  is  considered  to  be  a  deterministic  function  of  s,  0,  and  Therefore,  the  conditional 
pdf  of  the  data  is  related  to  the  pdf  of  n  by 

p(y|©,*,s)  =p„(y -sa(0,^)) .  (3) 

The  pdf  p(y  |0)  is  obtained  by  removing  variables  from  (3)  using  the  rule  for  conditional  pdf’s  [4] 


p(y|0)  =  j  j Pn(y-sa(0,^))p(s,'4'|0)  ds  d^. 

^  s 


(4) 


Substituting  (4)  into  (2)  yields  the  a  posteriori  pdf  of  0 

p(0|y)  =  J  j  p^{j-  sa(0,  ’*'))  p(s,  '®'|0)  ds  d^. 


(5) 


S 


We  can  simplify  (5)  by  applying  the  previous  assumptions  and  further  assuming  that  s,  and  0  are 
all  independent  of  each  other.  This  results  in  [3] 


p(0|y)  =  C'(y)p(0)  J  — ^J-^P(^) 


(6) 


where  C{y)  is  a  normalization  constant,  p(0)  is  assumed  uniform  over  the  parameter  space  of  interest,  and 

i?(0,'®')  =  |Fa"(0,^)R^iy|',  .. 

G(0,  =  Fa^^ (0,  ia(0,  '^)  +  1,  ^  ^ 

where  F  is  a  constant  and  Rn  =  F{nn^}.  The  maximum  a  posteriori  (MAP)  estimate  of  0  is  obtained  by 
maximizing  (6)  over  0,  i.e., 


©  =  argnmxp(©ly)  =  argn^ 


/ 


exn|HaSi\ 

G(0,'i') 


p(^)  d^. 


(8) 


If  the  vectors  a(0,^)  are  normalized  to  have  unit  norm  and  Rn  equals  the  identity  matrix,  equation  (8) 
can  be  written  as 

r  f  F^  Ia^f0  ’^)vl^l 

©  =  argn^  /  exp  < - - — ^  - —  >  p(^)  (9) 

For  numerical  implementation  of  this  estimator,  we  assume  that  p(’^)  is  a  uniformly  distributed  pdf  and 
approximate  the  integral  by  a  summation 


M 

0  =  arg  max  exp 

i=i 


F^|a"(0,»,)y|^ 

F+1 


I' 


(10) 


where  the  are  vectors  of  the  nuisance  parameters  sampled  from  their  probability  distributions.  A 
Monte  Carlo  approach  to  computing  (10)  was  proposed  in  [5],  Thus  M  inner  products  between  the  vectors 
,  M  and  y  are  computed  for  each  point  in  parameter  space  0. 
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I  Computationally-EfRcient  Maximum  A  Posteriori  Parameter  Es¬ 
timation 


The  exponential  in  (10)  can  be  approximated  over  a  finite  interval  by  a  linear  approximation  of  the  form 


exp  - 


F  +  1 


F  +  1 


+  b. 


(11) 


for  some  constants  a  >  0  and  b.  Substituting  this  approximation  into  (10)  gives 


aF^ 


M 


0  =  argmax  Mb  +  ^  |a^(©,  ^j)y\ 

j=i 


or  equivalently 


Af 


0  =  argmax^  la^(0,^^)y|' 


(12) 


(13) 


j-i 


Notice  that  if  the  perturbations  of  the  a(0,  vectors  for  each  0  over  the  realizations  are  small,  then 
the  variability  of  the  values  of  |a^(0,  will  also  be  small.  Under  this  condition,  (13)  will  be  a  close 

approximation  to  (10).  Rewriting  (13)  as 


0  =  axgii^y^Ra(0)y, 


(14) 


where 


M 


Ra(0)  =  2]a(0,'®'^)a^(0,’^,), 


(15) 


3^1 


reveals  that  the  assumption  of  small  a(0,  ^j)  vector  perturbations  over  implies  that  the  matrix  Ra(0) 
is  approximately  low-rank.  Therefore,  we  can  make  a  low-rank  approximation  to  Ra(€))  by  using  its  eigen 
decomposition  Ra(0)  =  USU^  and  setting  the  small  eigenvalues  in  S  to  zero  [6].  Defining  Si  as  a 
diagonal  matrix  of  the  non-zero  eigenvalue^and  Ui  as  a  matrix  of  the  corresponding  eigenvectors,  a  low- 
rank  approximation  to  Ra(0)  is  given  by  Ra(©)  =  UiSiUf^.  Substituting  this  approximation  into  (14) 
and  rewriting  the  expression  as  the  norm-squared  of  a  matrix-vector  product  gives 


0  =  axgn^  ||B; 


.(0)y||^ 


(16) 


where  Ba(0)  =  Sf  Uf . 

The  advantage  to  using  the  approximate  MAP  estimator  of  (16)  over  the  MAP  estimator  of  (10)  is  that 
the  computation  of  Ba(0)  can  be  done  off-line  before  the  processing  of  data.  The  number  of  on-line  inner 
product  computations  over  each  point  in  the  desired  parameter  space  0  is  reduced  and  equal  to  the  rank 
of  Ba(0).  In  contrast,  (10)  must  compute  M  inner  products  over  each  point  ©  on-line.  Thus,  (16)  is  a 
computationally-efficient  approximation  to  the  MAP  estimator  and  will  be  called  the  fast,  approximate  MAP 
estimator  or  FASTMAP. 


II  Simulation  Example  from  Matched-Field  Processing 

In  this  section  we  will  demonstrate  the  utility  of  the  FASTMAP  estimator  using  the  array  processing  tech¬ 
nique  known  as  matched-field  processing  (MFP)  [7].  MFP  is  a  model-based  source  localization  method  which 
is  a  function  of  the  ocean  environmental  parameters.  Precise  knowledge  of  the  environmental  parameters 


117 


are  required.  When  the  environmental  parameters  are  uncertain,  a  MAP  approach  to  MFP  can  be  used. 
Using  the  model  of  (1)  for  this  case,  ©  would  contain  the  source  location  parameters  of  interest  and  the 
uncertain  environmental  parameters  would  be  contained  in  The  vector  a(©,  is  called  a  replica  vector 
in  the  matched-field  literature. 

We  will  now  compare  the  outputs  of  the  MAP  estimator  and  the  FASTMAP  estimator,  through  a 
simple  simulation,  for  increasing  levels  of  environmental  uncertainty.  As  the  environmental  uncertainty 
level  increases,  the  perturbations  of  the  replica  vectors  at  each  0  over  environment  also  increase.  This  will 
illustrate  the  performance  decrease  of  the  FASTMAP  estimator  when  the  small  perturbation  assumption  is 
no  longer  valid. 

We  simulated  a  typical  shallow  water  environment.  All  of  the  environmental  parameters  were  assumed 
known  except  one.  A  linear  sound  velocity  profile  was  used  which  had  a  sound  speed  at  the  surface  of  1500 
m/s  and  a  sound  speed  at  the  bottom  of  1480  m/s.  We  introduced  uncertainty  by  allowing  the  sound  speed 
at  the  bottom  to  be  uncertain.  A  20-element  vertical  array,  which  spanned  the  water  column,  was  used. 
The  data  vector  y  was  generated  using  the  nominal  value  of  the  sound  speed  at  the  bottom  1480  m/s,  while 
the  estimators  assumed  its  value  could  lie  anywhere  in  its  specified  uncertainty  interval  of  1480  ±  As  m/s. 
The  source  was  located  at  a  range  of  9000  m  and  a  depth  of  24  m.  For  convenience,  it  was  assumed  that 
the  source  depth  was  known  a  priori  by  both  processors.  Therefore,  range  was  the  only  source  location 
parameter  each  had  to  search  over.  Figures  1  and  2  show  the  outputs  of  both  estimators  for  As  =  1, 7.5 
m/s,  respectively.  For  these  levels  of  uncertainty,  both  estimators  are  able  to  correctly  estimate  the  source 
range.  Figure  3  shows  the  outputs  as  As  was  increased  to  25  m/s.  At  this  point,  the  small  perturbation 
assumption  no  longer  holds.  With  this  large  degree  of  uncertainty,  the  FASTMAP  estimator  produces  a 
peak  at  an  incorrect  range  estimate,  while  the  MAP  estimator  displays  a  peak  at  the  correct  estimate. 

This  simple  simulation  provides  an  example  of  the  computational  savings  of  using  the  FASTMAP  esti¬ 
mator  over  the  MAP  estimator.  W^ith  As  =  7.5  m/s,  the  uncertainty  interval  of  the  sound  speed  at  the 
bottom  ranged  from  1472.5  m/s  to  1487.5  m/s.  This  interval  was  sampled  at  0.5  m/s  resulting  in  30  sample 
points,  i.e.,  M  =  30.  However,  the  rank  of  Ba(0)  was  equal  to  5.  Therefore,  only  5  inner  products  per  0 
were  required  for  the  FASTMAP  estimator,  while  the  MAP  estimator  required  30. 

A  more  complex  MFP  example  in  which  seven  environmental  parameters  were  uncertain  is  shown  in  [8]. 
Using  a  smoothed  version  of  the  FASTMAP  estimator,  it  is  shown  that  the  FASTMAP  estimator  provides 
a  computational  savings  of  75:1  over  the  MAP  estimator. 


Ill  Conclusion 

The  maximum  a  posteriori  (MAP)  parameter  estimator  was  derived  for  a  data  model  applicable  to  array 
processing  and  time-series  analysis.  It  was  then  shown  that  a  computationally-efficient  approximate  MAP 
estimator,  called  the  FASTMAP  estimator,  could  be  obtained  if  the  perturbations  of  the  model  vectors 
a(©,  ^j)  at  each  point  in  parameter  space  0  over  the  nuisance  parameter  space  were  small.  Simulation 
results  from  matched-field  processing  were  shown  to  illustrate  the  effects  on  the  FASTMAP  estimator  when 
the  small  perturbation  assumptions  are  violated. 
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ts  of  FASTMAP  and  MAP  processors  in  the  presence  of  2  m/s  sound-speed  uncertainty. 
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Outputs  of  FASTMAP  and  MAP  processors  in  the  presence  of  15  m/s  sound-speed  uncertainty. 
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Figure  3:  Outputs  of  FASTMAP  and  MAP  processors  in  the  presence  of  50  m/s  sound-speed  uncertainty. 
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ABSTRACT 

Multiple-frequency  robust  ad^tive  matched-field  processing  (MFP)  has  been  applied  to  data 
from  both  the  Hudson  Canyon  Experiment'*^  and  the  Elba  Sea  Experiment®.  The  Hudson 
Canyon  data  set  consisted  of  20  snapshots  of  multiple-cw  transmissions  and  was  distributed  in 
the  5th  MFP  workshq>.  The  Elba  Sea  data  set  consisted  of  a  period  of  broadband  M-sequence 
transmissions  and  was  distributed  in  the  6th  MFP  workdiop.  The  MFP  results  from  both  data 
sets  demonstrate  excellent  source  localization  up  to  about  5  km  in  range  in  shallow  water.  In 
a  November  1995  experiment  in  the  Gulf  of  Mexico,  a  vertical  array  was  deployed  in  water  188 
m  deep  .  During  a  TL  run,  a  projector  transmitting  multiple-cw  tones  was  towed  along  the 
isobath  line  between  5.5  km  and  15.5  km  away  from  the  array.  The  MFP  results  show  that 
multiple-frequency  robust  adaptive  processing  provides  excellent  source  localization  up  to  15  km 
in  this  shallow  water  environment. 


1.  INTRODUCTION 

MFP  refers  to  array  processing  algorithms  which  exploit  the  full  field  structure  of  signals 
prq)agating  in  an  ocean  waveguide.  In  it  an  acoustic  model  is  used  to  calculate  replica  vectors 
in  a  search  space.  The  measured  signal  plus  noise  data  vectors  are  then  matched  wito  the  replica 
vectors  and  an  ambiguity  surface  is  computed.  The  locations  of  peaks  in  the  ambiguity  surface 
are  estimates  of  the  source  location.  An  excellent  overview  of  MFP  can  be  found  in  the  paper 
by  Baggeroer,  et.  al. 

Adaptive  matched  field  processing  (MFP)  is  an  q}timal  method  for  source  localization  and 
detection  in  the  ocean.  It  uses  the  measured  signal-plus-noise  (data)  vectors  to  minimize  the 
contributions  from  those  components  that  do  not  mamh  with  the  steering  vector  to  a  given  search 
cell.  When  there  is  a  mismatch  between  the  actual  signal  vector  and  the  calculated  or  replica 
vector,  there  is  a  degradation  in  the  performance  of  the  adaptive  processor  which  can  be  severe 
when  the  signal-to-noise  ratio  (SNR)  is  high. 

A  review  of  several  traditional  approaches  as  well  as  some  newer  algorithms  to  make 
adaptive  MFP  robust  to  mismatch  has  been  discussed  by  Lee,  et.al.®.  Applications  of  robust 
adaptive  algorithms  to  multiple-frequency  signals  has  been  shown  by  Lee®  using  the  data  set 
from  the  Hudson  Canyon  Experiment.  Comparing  with  arithmetic  mean  (averaging  the 
amplitudes),  taking  the  geometric-mean  (averaging  dB's)  or  the  harmonic-mean  (averaging  the 
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reciprocal  of  the  amplitudes)  of  the  discrete  frequency  output  of  the  robust  adaptive  algorithms 
for  multiple-frequency  signals  have  demonstrated  excellent  source  localization  results. 

In  this  study,  two  robust  adaptive  algorithms,  namely  the  Feedback-Loop  White-Noise- 
Constrained  (FLWNC)  method  and  the  Signal-Coherence-Constrained  Reduce-Rank  (SCCRR) 
method  have  been  ^plied.  Section  2  will  review  both  FLWNC  and  SCCRR  algorithms.  In 
section  3  we  will  describe  the  setup  of  a  recent  shallow  water  experiment  off  Gulf  of  Mexico, 
in  section  4  we  will  show  MFP  results,  and  in  section  5  we  will  present  a  discussion  of  the 
results. 

2.  MULTIPLE-FREQUENCY  ROBUST  ADAPTIVE  ALGORITHMS 

Adaptive  processing  uses  the  measured  signal  plus  noise  data  vectors  to  minimize  the 
sidelobe  contributions  from  those  components  that  do  not  match  with  the  replica  vector  for  a 
given  search  cell.  Most  adjq)tive  techniques  used  today  had  their  genesis  in  the  Minimum 
Variance  Distortionless  Filter  (MVDF)  algorithm,  or  Capon's  Maximum  Likelihood  Method 
(MLM).  Let  K  be  the  measured  covariance  matrix;  from  singular  value  decomposition,  it  can 
be  decomposed  into  a  set  of  eigenvectors  V  associated  with  a  set  of  eigenvalues  Xi  .  Let  A 
be  a  replica  vector  for  a  given  search  cell  (look  direction).  The  MVDF  minimizes  the  variance 
at  the  output  of  a  linear  weighing  of  the  sensors  subject  to  the  distortionless  constraint  that 
signals  in  the  "look  direction"  have  unity  gain.  Its  ouq)ut  is  given  by 


Without  mismatch,  the  steering  vector  A  is  perfectly  matched  with  signal  eigenvectors  in  the 
"look  direction",  the  rest  of  eigenvectors  are  orthogonal  to  the  steering  vector  and  have  no  effect 
on  the  signal  output.  But,  when  there  is  mismatch  present,  the  steering  vector  is  no  longer 
orthogonal  to  the  rest  of  eigenvectors.  The  noise  vectors  associated  with  the  least  signifrcant 
eigenvalues  then  dominate  the  inverse  processing  in  Eq.(l)  and  degrade  the  signal  estimation. 

In  MFP,  there  are  several  forms  of  mismatch,  including  mismatch  due  to  environmental 
uncertainties  and  fluctuation  and  mismatch  due  to  system  errors.  All  forms  of  mismatch  are 
either  deterministic  or  random.  Deterministic  mismatch  degrades  the  signal  estimation  and 
causes  localization  bias,  but  it  can  be  minimized  if  more  ground  truth  is  provided.  Random 
mismatch  degrades  the  signal  estimation  but  will  not  bias  the  localization.  Random  mismatch 
cannot  be  minimized  so  that  robust  algorithms  were  developed  to  tolerate  a  certain  level  of 
random  mismatch. 

The  FLWNC  method  reduces  sensitivity  to  mismatch  by  adding  white  noise  (e)  to  the 
diagonal  elements  of  the  covariance  matrix  K.  Adding  white  noise  (e)  to  the  diagonal  elements 
of  the  covariance  matrix  K  is  the  same  as  adding  e  to  each  eigenvalue  without  modifying  the 
eigenvectors.  The  white-noise-constrained  output  that  results  is 

For  each  search  location,  the  MVDF  white-noise  processing  gain  and  spectral  output  are 
calculated.  Tfre  white-noise  processing  gain,  defined  as  the  amplitude  squared  of  the  weight 
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(2) 


'me 


\vljf 


vector,  is  directly  proportional  to  the  SNR.  If  the  white-noise  processing  gain  falls  below  a 
constraining  value,  a  "noise-cell"  is  defined  and  the  FLWNC  spectral  oulput  is  set  equal  to  the 
MVDF  spectral  output.  If  the  white-noise  processing  gain  is  above  the  constraining  value,  a 
"signal-cell"  is  defined  and  an  amount  of  white  noise  that  equals  the  MVDF  q)ectral  output  is 
added  in  the  processing  and  a  new  white-noise  processing  gain  and  a  new  q)ectral  output  are 
calculated.  TTie  process  is  repeated  until  the  white-noise  processing  gain  falls  below  the 
constraining  value.  The  feedback-loop  ^proach  insures  that  the  amount  of  white  noise  added 
in  each  iteration  is  adequate  so  that  the  iteration  procedure  converges  rapidly  without 
overshooting. 

For  each  search  location,  the  coherence  between  an  eigenvector  and  the  steering  vector, 
defined  as  the  dot  product  of  these  two  unit  vectors,  is  calculated.  Consider  an  environment  in 
which  the  signal  coherence  is  bounded  within  a  constraining  value  a.  In  the  SCCRR  method, 
the  coherence  is  compared  with  the  constraining  value  a.  If  it  is  greater  than  a  the  eigenvector 
falls  in  a  signal  ^ace  for  this  search  location  otherwise  the  eigenvector  falls  in  a  noise  space. 
The  SCCRR  method  excludes  eigenvectors  in  the  noise  space  and  uses  those  in  signal  space  for 
the  processing.  The  SCCRR  ouq>ut  is  then 


Broadband  MFP  is  calculated  in  general  by  taking  the  weighted  or  unweighted  arithmetic- 
mean  of  ambiguity  surfaces  resulting  from  each  narrowband  component  across  the  frequency 
band  of  interest.  Instead  of  taking  the  arithmetic  mean,  taking  the  geometric-mean  or  the 
harmonic-mean  has  proven  to  be  a  more  effective  way  to  control  the  ambiguous  sidelobes.  Let 
y(fi )  be  the  narrowband  MFP  output  at  frequency  f  ,  and  ¥»,  Y,  ,Yi  be  the  arithmetic-mean, 
geometric-mean,  and  harmonic-mean  over  frequency,  req)ectively.  Mathematically,  Y»  >  Yj 
>  Yh ;  they  become  equal  when  y(fi  )  are  the  same.  The  narrowband  responses  line  up  at  the 
target  location  and  all  three  means  have  the  same  result.  Because  sidelobe  distributions  are 
different  (misaligned)  among  the  narrowband  re^onses,  the  geometric-mean  and  the  harmonic- 
mean  have  lower  sidelobes  than  the  arithmetic-mean. 

3.  EXPERIMENTAL  SETUP 

A  shallow  water  experiment  was  conducted  in  November  1995  in  Gulf  of  Mexico 
approximately  100  miles  directly  south  of  Panama  City.  A  vertical  array  was  deployed  at 
28”23.063  N  SS^IT-SST  W.  The  water  depth  was  measured  to  be  approximately  188  m  at  the 
array  site.  The  vertical  array  had  32  hydrophones  and  a  total  aperture  of  153.  The  upper  13 
hydrq)hones  were  uniformly  i^aced  at  8  m,  and  the  lower  19  hydrophones  were  uniformly 
spaced  at  3  m.  It  was  deployed  in  a  bottom  moored  configuration  and  the  deepest  hydrophone 
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was  3  m  above  the  bottom.  Few  hydrophones  were  not  functional  and  were  excluded  in  MFP 
processing.  Although  array  navigation  data  were  collected  during  the  experiment,  they  were  not 
available  at  the  time  of  data  processing.  So  the  array  was  assumed  to  be  perfectly  vertical  for 
all  MFP  processing. 

During  the  experiment,  a  HLF-6  source  was  towed  at  a  speed  of  ^proximately  5  knots  for 
a  transmission  loss  measurement.  The  HLF-6  source  brmdcasted  5  CW  tones  at  47  Hz,  97  Hz, 
147  Hz,  197  Hz,  and  247  Hz  at  a  nominal  source  level  of  195  dB  per  line.  Figure  1  shows  the 
track  where  the  source  was  towed  starting  from  deep  to  shallow  and  then  towed  along  an  isobath 
line.  In  this  study  we  concentrate  on  the  source  localization  along  the  isobath  line. 

4.  EXPERIMENTAL  RESULTS 

Figure  2  shows  the  average  XBT  profile  taken  close  in  time  for  this  measurement.  The 
sediment  properties  were  modified  from  the  Nobles  parameters®  to  provide  the  best  results  for 
all  frequencies  at  all  ranges.  A  range-independent  normal  mode  code  was  used  for  matched-field 
processing.  Only  the  97  Hz,  147  Hz,  and  197  Hz  signals  were  processed  and  the  geometric- 
mean  results  of  tiiese  three  frequencies  are  presented. 

Figure  3  diows  the  range-depth  responses  of  the  conventional  MFP  for  the  source  at  ten 
different  ranges  from  5.63  km  to  14.45  km.  Each  range-depth  response  contour  is  normalized 
to  its  peak  and  has  a  6  dB  dynamic  range.  Sidelobes  tqipear  at  about  3  to  4  dB  below  the  peak. 
Figure  4  shows  the  range-depth  responses  of  the  FLWNC  processor  and  Figure  5  shows  the 
range-depth  responses  of  the  SCCRR  processor.  A  40  dB  dynamic  range  is  used  for  these 
response  contours.  Sidelobes  appear  at  about  20  dB  below  the  peak  for  the  FLWNC  processor 
and  at  more  than  30  dB  below  the  peak  for  the  SCCRR  processor.  These  results  show  excellent 
source  localization  at  all  ranges. 

Figure  6  shows  signal  mismah;h.  Signal  mismatch  is  defined  as  the  difference  between  the 
array  received  energy  and  the  peak  ouq)ut  in  the  conventional  MFP.  Along  the  isobath  line, 
from  5  km  to  15  km,  the  signal  mismatch  is  almost  independent  of  range  and  is  about  1  dB  for 
97  Hz,  1.3  dB  for  147  Hz,  and  2  dB  for  197  Hz.  When  &e  source  moved  from  deep  to  shallow 
water,  at  ranges  from  20  km  to  15  km,  signal  mismatch  increased.  Figure  7  shows  the 
localization  performance  for  different  frequencies  at  different  ranges.  The  source  localization  is 
consistent  for  all  three  frequencies  along  the  isobath  line  up  to  15  km.  The  consistency  degrades 
when  the  distance  is  greater  than  15  km  because  the  water  depth  changed  which  is  not  included 
in  our  range-independent  MFP  processing.  A  consistent  bias  in  depth  estimation  indicates  that 
the  source  may  not  have  been  towed  at  the  planned  depth. 

5.  CONCLUSIONS 

Multiple-frequeiKjy  robust  adaptive  matched-field  processing  has  been  applied  to  data 
collected  from  a  shallow  water  test  in  Gulf  of  Mexico.  Using  the  average  profile  collected  close 
in  time  to  the  measurement,  the  range-independent  MFP  provides  an  excellent  source  localization 
up  to  15  km  along  an  isobath  line. 
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ABSTRACT 

The  goals  of  this  research  is  to  use  MFP  technique  to  map  the  spatial  distribution  of 
reverberation  as  a  tool  for  understand  the  physics  of  ocean  surface  scattering.  In  a  November 
1995  experiment  in  the  Gulf  of  Mexico,  a  vertical  array  was  deployed  in  water  188  m  deep  and 
a  CW  source  was  placed  on  the  bottom  760  m  away  from  the  array.  The  received  energy  was 
qiectrally  separated  from  the  carrier  to  isolate  the  components  scattered  from  the  sea  surface  and 
then  MFP  was  performed  to  examine  its  spatial  distribution.  A  simulation  was  performed  by 
generating  time  step  realizations  of  an  ocean  surface  with  randomly  distributed  surface  waves. 
A  modified  FEPE  to  include  rough  surface  was  used  to  propagate  signals  to  a  vertical  array.  In 
addition  to  the  surface  forward-scattered  energy  which  is  modeled  in  the  simulation,  bottom- 
back-scattered  energy  is  also  observed  in  the  MFP  data  reqionses. 

1.  INTRODUCTION 

MFP  combines  an  acoustic  propagation  model  with  matched-filter  processing  to  localize  a 
source  in  three-dimensions.  It  has  been  shown  to  be  effective  in  discriminating  passive  sources 
both  in  range  and  depth  by  exploiting  the  vertical  structure  of  the  signal.  In  active  sonar 
systems,  for  frequencies  below  1  kHz,  most  reverberation  comes  from  the  surface  and  bottom. 
MFP  can  discriminate  in  depth  and  should  have  potential  to  provide  significant  gain  against 
active  reverberation.  In  a  numerical  study,  Dozier,  et.  al.^‘^  showed  that  it  has  the  potential  to 
yield  gain  against  reverberation  not  far  from  the  ideal  array  gain.  We  have  applied  MFP 
technique  to  map  the  forward-scattered  signal  around  the  active  direct  blast®.  MFP  successfully 
maps  the  scattered  energy  near  the  surface  in  both  real  data  and  simulations. 

In  general,  the  active  scattering  problem  should  include  two  prq)agation  paths:  source  to 
scatterer  and  scatterer  to  receiver.  A  scattering  kernel  is  included  to  descril^  the  scattering 
mechanism.  In  a  stratified  medium,  the  acoustic  field  scattered  from  a  scatterer  can  be  modeled 
as 
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where  are  the  mode  functions  and  k„  are  their  associated  wave  numbers.  S  is  the  scattering 
kernel  that  describes  how  energy  in  one  mode  is  re-distributed  to  other  modes.  In  our  study,  we 
assume  that  the  scatterer  is  illuminated  by  the  source  and  "re-radiates"  energy  isotropically.  In 
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The  prq)agation  path  from  source  to  scatterer  affects  only  the  amplitude,  but  not  the  structure 
of  the  scattered  signal.  Hie  scattered  acoustic  field  received  at  a  vertical  array  depends  only  on 
the  propagation  path  from  scatter  to  receiver.  We  then  could  use  passive  matched-field 
techniques  to  localize  all  the  scatterers  and  characterize  their  ^atial  distribution. 

To  avoid  the  complication  of  long-range  multiple  bounces,  in  this  paper  we  study  MFP 
responses  of  one-  and  two-bounce  scattered  signals.  The  goal  is  to  characterize  the  ^atial 
distribution  of  the  scattered  field  and  to  better  understand  the  scattering  mechanism  in  the 
temporally  and  ^atially  varying  ocean  environment. 

In  section  2,  we  will  describe  the  setup  of  a  1995  shallow  water  experiment.  In  section  3, 
a  simulation  using  a  rough  surface  FEPE  will  be  described.  Results  of  real  data  and  simulation 
will  be  presented  in  section  4  and  a  summary  will  be  included  in  section  4. 

2.  EXPERIMENTAL  SETUP 

In  November  1995,  a  test  was  conducted  in  the  northeastern  Gulf  of  Mexico.  The  objective 
of  the  test  is  to  collect  data  for  validating  models  for  predicting  reverberation  distributions.  One 
of  the  concqpt  to  be  validated  is  how  well  MFP  techniques  can  be  used  to  study  reverberation 
distribution  in  range  and  depth.  The  test  site  and  season  had  been  selected  to  minimize  the 
complexity  of  2u;oustic  modeling  while  still  providing  a  wide  range  of  sea  conditions. 

A  vertical  array  was  deployed  in  water  188  m  deep.  The  array  consisted  of  32  hydrqihones 
with  the  tq)  hydrqihone  at  a  depth  of  32  m;  the  upper  13  hydrqihones  were  uniformly  ^aced 
at  8  m,  and  the  lower  19  hydrqihones  were  uniformly  spaced  at  3  m.  Although  the  array  was 
monitored  by  a  tran^nder  system,  no  array  element  location  data  are  available  at  this  time. 
So,  the  array  was  assumed  to  be  perfectly  vertical  for  this  study.  A  CW  source  was  placed  on 
the  bottom  a  short  distance  away  from  the  vertical  array.  It  broadcasted  a  240.375  Hz  CW  with 
a  source  level  of  qiproximately  175  dB.  Hie  CW  source  was  localized  acoustically  using  MFP 
at  760  m  from  the  vertical  receiving  array. 

Hie  plan  was  to  turn  the  source  on  for  2  to  3  hours  continuously  4  times  a  day  for  5 
consecutive  days  to  sample  various  sea  conditions.  But  due  to  the  difficulty  experienced  during 
the  experiment,  only  10  minutes  data  were  acquired  for  the  short  range  measurement.  Hie  wind 
speed  was  about  10  knots  and  the  significant  wave  height  was  about  2  feet.  More  than  twenty 
sound  speed  profiles  were  collected,  using  XBTs  and  CTDs,  within  the  two  days  proceeding  the 
acoustic  data  collected.  Figure  1  shows  the  average  sound  qieed  profile  which  is  used  in  MFP 
processing.  The  sediment  properties  were  modified  from  Knobles  parameters®  to  best  localize 
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a  towed  source  at  various  ranges. 

3.  SIMULATION 

The  finite-element  parabolic  equation  (FEPE)  method  for  acoustic  propagation  has  been 
extended  to  handle  a  rough  air-water  interface  treated  as  a  series  of  stair  steps^^\  Because  finer 
vertical  spacing  is  needed  to  resolve  the  interface,  finer  vertical  grids  near  the  surface  and  coarser 
grids  in  die  rest  of  the  ocean  are  used.  In  our  simulation,  below  3  times  of  the  maximum  surface 
wave  height  a  coarse  spacing  of  2.5  m  is  used,  and  above  a  finer  spacing  of  8  times  denser  is 
used. 

A  rough  sea  surface  moving  parallel  to  the  source-receiver  direction  with  a  (20-knot) 
Pierson-Moskowitz  spectrum  was  projected  onto  the  surface  of  the  2D  vertical  plane  containing 
the  source  and  receiver.  Simulation  was  done  by  generating  correlated  realizations,  every  half 
seconds,  of  an  ocean  surface  with  randomly  distributed  surface  waves  and  using  the  modified 
rough  surface  FEPE  to  prqiagate  signals  to  the  array. 

4.  RESULTS 

Phone  spectra  from  data  and  simulation  were  calculated  using  an  FFT  with  Hanning  shading. 
Figure  2  shows  the  array  average  spectrum  from  data.  We  separated  the  Doppler  spectrum  into 
three  bands,  the  Carrier,  Down-Doppler,  and  Up-Dqipler,  and  applied  MFP  using  a  range- 
independent  normal  mode  code  to  each  band.  Figure  3  shows  the  range-depdi  responses  of  the 
Carrier  frequency.  The  experimental  results  are  shown  in  the  top  two  panels  and  the 
corresponding  simulation  results  are  shown  in  the  bottom  two  panels.  The  conventional  MFP 
re^onses  ate  shown  on  the  left  and  the  robust  adaptive  MFP  responses^*^  are  shown  on  the  right. 
TTie  source  was  localized  at  760  m  on  the  bottom.  The  conventional  MFP  responses  show  high 
ambiguous  sidelobes  which  depict  the  propagation  paths  outward  from  the  source.  It  shows  two 
paths  reach  the  array:  one  path  direct  propagates  to  the  array  and  the  other  path  hits  the  surface 
and  is  refracted  to  die  array.  The  robust  adaptive  MFP  suppresses  the  sidelobes  and  shows  an 
excellent  localization  of  the  source.  Figure  4  shows  the  MFP  rehouses  of  the  down-Doppler 
energy  and  the  Up-Dc^pler  energy.  The  forward-scattered  energy  is  localized  near  the  place 
where  the  acoustic  prqiagation  paths  reach  the  surface.  The  rough-surface  FEPE  simulated 
spatial  distributions  of  die  forward-scattered  energy  are  consistent  with  the  real  data  results.  In 
additional  to  the  forward-scattered  energy  other  scattered-components  which  might  be  related  to 
the  bottom-back-scattered  energy  and  volume-scattered-energy  are  observed  in  the  MFP  responses 
of  the  real  data. 

5.  SUMMARY 

In  this  study  we  have  demonstrated  that  MFP  technique  can  be  used  to  study  the  range-depth 
spatial  distribution  of  reverberation.  More  information  can  be  obtained  to  better  understand  the 
physics  of  how  the  ocean  surface  is  illuminated  and  re-radiates  the  energy. 
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Abstract 


The  effectiveness  of  novel  regularization  techniques  for  generating  stable  solutions  to  the  ill-conditioned 
linear  problem  of  modal  decomposition  was  examined.  These  techniques  were  based  on  either  a  specific 
Tikhonov  regularization  scheme  or  one  of  several  more  general  spectral  regularizing  algorithms.  The 
model  system  for  this  simulation  study  consisted  of  an  ideal  waveguide  and  a  vertical  array  spanning  half 
the  water  column.  It  was  verified  that  regularization  was  necessary  to  obtain  a  stable  solution,  and  that 
the  accuracy  of  estimation  depended  on  the  regularization  scheme  and  the  nature  of  the  noise.  When 
the  noise  was  assumed  to  be  known  exactly  (i.e.,  the  deterministic  case),  spectral  regularization  gave 
the  most  accurate  result.  For  the  situation  where  only  statistical  knowledge  of  the  noise  was  available, 
Tikhonov’s  scheme  was  somewhat  more  accurate  than  the  spectral  approach.  The  accuracy  of  the  reg¬ 
ularized  solution  was  observed  to  depend  strongly  on  the  characteristics  of  the  right-hand  side  of  the 
system  of  linear  equations  as  well  as  on  the  condition  of  the  matrix.  The  most  accurate  solutions  were 
obtained  when  the  source  depth  was  within  the  depth  region  spanned  by  the  array. 


I  Modal  Decomposition 

Modal  decomposition  is  an  essential  stage  in  the  use  of  matched-mode  processing  (MMP)  techniques  for 
source  detection  and  localization.  MMP  can  be  used  in  conjunction  with  a  nearest  neighbours  algorithm  to 
yield  processing  techniques  which  are  substantially  more  efficient  than  matched-field  processing.  MMP  first 
requires  the  estimation  of  the  m  complex  mode  amplitudes  x  from  the  measured  fields  y  at  the  n  sensors, 
by  solving  the  n  x  m  system  of  linear  equations  Ax  =  y,  where  A  is  the  design  matrix.  We  shall  consider 
here  only  the  case  of  vertical  arrays,  for  which 
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where  (l)j{zk)  is  the  mode  function  amplitude  for  the  jth  mode  and  the  fcth  sensor  and  x  is  a  vector  of  mode 
excitations  at  the  source  depth  i.e., 
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where  (*)  denotes  the  conjugate  transpose. 


II  Inversion  and  Regularization 

Given  A  and  y,  where  the  elements  of  y  are  the  noise- containing  measured  field  values,  the  least-squares 
estimate  of  the  mode  amplitudes  x  is  characterized  by  the  normal  equations: 

A*  Ax  =  A’^y. 

Provided  a  unique  solution  exists,  it  can  be  conveniently  expressed  in  terms  of  the  eigendecomposition  of 
A’^A;  i.e., 


EDE*x  =  A*y,  from  which 

X  =  ED“^E*A*y, 


where 

E  is  a  matrix  of  the  eigenvectors  of  A*^  A, 

D  is  a  diagonal  matrix  with  diagonal  elements  A* , 

D-1  is  a  diagonal  matrix  with  diagonal  elements  l/Ajt,  and 
Ajk  is  the  kth  eigenvalue  of  A*  A. 

For  the  case  where  the  array  spans  only  a  portion  of  the  water  column,  the  matrix  A*  A  can  be  very 
ill-conditioned.  For  example,  for  the  case  of  an  array  spanning  only  the  top  half  of  a  100-m  deep  waveguide, 
the  condition  number  of  this  matrix  was  found  to  increase  from  10^  to  5  x  10^^  as  the  number  of  modes 
was  increased  from  5  to  10.  This  resulted  in  unstable  and  highly  inaccurate  estimates  of  x  when  the  normal 
equations  solution  was  computed. 

The  inversion  (i.e.,  estimation  of  x)  may  be  stabilized  through  regularization.  Regularizing  operators  Rs 
provide  a  stabilized  approximation  to  the  inverse  operator,  and  ensure  that  solutions  obtained  by  inversion 
converge  to  the  true  solutions  as  errors  in  the  measurements  tend  to  zero.  Design  of  suitable  regularizing 
operators  depends  on  the  nature  of  the  prior  knowledge  about  the  signal  and  noise  characteristics.  These 
operators  can  be  expressed  in  spectral  form  in  terms  of  the  eigendecomposition  above,  where  dkkj  the  diagonal 
elements  of  the  matrix  (the  regularized  approximation  to  D”^),  are 

dkk  =  >^fc/(A|  +  ^^(Ajk)), 

where  ^^(Ajk)  is  a  nonnegative  stabilizing  function. 

Standard  Tikhonov  regularization  [1]  is  a  global  scheme  which  involves  minimizing 

l|Axa-y|l2  +  a||x||2, 
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for  which  the  form  of  the  approximating  operator  is 


=  (A*A  +  aI)-^A*.  (1) 

We  consider  here  the  specific  Tikhonov  regularizing  operator  which  is  obtained  by  finding  a  —  a (6)  and 
to  satisfy  the  discrepancy  equation: 

G{a)  =  ||Axa  -  y|p  -  =  0, 

where  is  a  constant  which  depends  on  the  noise  level  and  model  errors, 

Spectr^d  Regularization  (SR)  [2]  is  a  more  flexible  form  of  regularization  than  the  Tikhonov  scheme. 
It  may  be  expressed  in  the  following  form: 

=  Y^WkGk, 

k 

where 

Cjb  is  the  ikth  eigenvector  of  A*  A, 

uj^  =  (e^  '  A* y)  is  the  A:th  component  of  the  spectral  decomposition  of  A'^y,  and 
Wk  is  the  kth  spectral  component  of  the  solution  x. 

We  note  that  for  the  standard  Tikhonov  regularizing  operator  (Eq.  1),  qsi^k)  =  Thus  this  operator 
is  a  simple  form  of  spectral  regularizing  operator.  However,  in  the  following  we  shall  restrict  the  use  of  the 
term  “spectral  regularization”  to  those  forms  described  below. 

In  determining  appropriate  forms  for  qs  and  w,  one  approach  is  to  assume  that 

\uk  -  Wil  <  0-k- 

In  this  case  it  may  be  shown  that  the  solution  which  minimizes  |u;jfe|  is 


P  J  0  if  iui  I  <  cTi 

“'""I  • 

This  form  is  termed  the  initial  solution. 

Once  this  initial  solution  has  beed  determined,  the  following  expression  may  be  iterated: 

,,,n+l  _  ^kUk 


On  convergence  (i.e.,  when  n  — »  oo),  this  form  is  termed  the  limit  solution. 

Both  forms  require  that  an  estimate  of  the  spectral  decomposition  of  the  noise  ak  be  provided  for  each 
spectral  component  k.  This  was  done  using  one  of  two  noise  models:  (1)  deterministic ,  where  the  known 
noise  values  n  at  the  sensors  were  used:  i.e.,  set 


<^k  =  |(e%  ■  A*n)|, 

and  (2)  statistical,  where  the  RMS  spectral  decomposition  of  the  noise  was  estimated  for  various  realizations 
of  noise  iij  and  scaled  by  a  constant  c:  i.e.,  set 

O-jfe  =  c  X 
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Here,  c  was  set  to  2.0  (twice  the  RMS  noise  level);  this  value  was  chosen  to  suppress  the  effects  of  the  larger 
random  realizations  of  noise  in  the  spectral  components  of  the  data.  Thus,  only  when  the  signal-hnoise  spec¬ 
tral  level  in  the  data  was  at  least  twice  the  estimated  average  noise  level  did  the  eigenvector  corresponding 
to  that  spectral  component  contribute  to  the  solution. 


Ill  Comparison  of  Regularization  Schemes 

The  regularization  schemes  described  above  were  compared  in  a  series  of  simulations  involving  an  ideal 
waveguide.  The  waveguide  was  100  m  deep,  with  a  pressure-release  surface  and  a  rigid  bottom.  A  vertical 
array  was  positioned  so  as  to  span  the  top  or  bottom  half  of  the  waveguide,  and  a  source  at  a  frequency 
between  15  and  75  Hz  (giving  between  2  and  10  modes)  was  modelled  at  a  specified  depth  and  range.  The 
accuracy  of  the  computed  solution  was  expressed  in  terms  of  the  relative  sum  of  squares  of  deviations:  i.e., 

misestimated  “  ^true  I  i  /Untrue 1 1  • 

The  following  observations  were  made  as  a  result  of  these  simulations: 

1.  The  SR  scheme  with  deterministic  noise  was  more  accurate  than  the  Tikhonov  scheme  (but  this  type 
of  prior  knowledge  about  the  noise  is  not  available  in  applications  such  as  the  present  one,  where  the 
noise  is  random).  The  SR  scheme  with  statistical  noise,  though  realizable  in  practice,  was  somewhat 
less  accurate  than  the  Tikhonov  scheme. 

2.  The  initial  form  for  SR  was  slightly  more  accurate  than  the  limit  form. 

3.  There  was  very  little  effect  on  accuracy  of  increasing  the  number  of  sensors  in  the  array. 

4.  The  accuracy  tended  to  decrease  rapidly  as  the  number  of  modes  was  increased,  but  this  decrease  was 
not  monotone. 

5.  The  accuracy  depended  on  range,  but  again  the  dependence  was  not  monotone. 

6.  The  accuracy  depended  on  the  depth  at  which  the  source  was  placed  in  the  channel. 

This  last  observation  is  illustrated  in  Figures  1  and  2,  where  the  relative  deviations  of  the  solution  were 
estimated  as  a  function  of  the  depth  of  the  source.  In  these  simulations,  the  frequency  was  40  Hz  (resulting 
in  five  modes)  and  the  source- array  range  was  2000  m.  Figure  1  shows  the  case  for  an  array  in  the  top  half 
of  the  water  column,  and  Figure  2  shows  the  results  for  an  array  in  the  lower  half.  It  is  evident  that  the 
most  accurate  solutions  were  obtained  under  conditions  where  the  source  depth  was  within  the  depth  region 
spanned  by  the  array. 

This  observation  is  also  of  theoretical  interest,  since  under  these  conditions  the  matrix  A  for  a  particular 
array  configuration  is  constant  as  the  source  depth  varies.  Therefore  the  accuracy  of  the  solution  is  strongly 
influenced  by  the  characteristics  of  the  right-hand  side  of  the  system  of  equations  as  well  as  by  the  condition 
number.  The  physical  explanation  underlying  this  mathematical  result  is  an  area  for  future  investigation.^ 
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Array  Elements;  0  —  50  m 


Figure  1:  Relative  sum-of-squares  (SSQ)  of  deviations  of  regularized  solutions  as  a  function  of  source  depth, 
for  an  array  in  the  top  half  of  the  channel.  These  results  were  obtained  for  a  40- Hz  source  in  a  100-m  channel 
(5  modes)  and  a  source-array  range  of  2000  m. 


Array  Elements:  50—100  m 


Source  Depth  (m) 

Figure  2:  Relative  sum-of-squares  (SSQ)  of  deviations  of  regularized  solutions  as  a  function  of  source  depth, 
for  an  array  in  the  lower  half  of  the  channel.  These  results  were  obtained  for  a  40-Hz  source  in  a  100-m 
channel  (5  modes)  and  a  source- array  range  of  2000  m. 
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Abstract 

This  paper  applies  the  method  of  matched-field  processing  (MFP)  to  localize  multiple 
sources  in  cases  where  the  bathymetry  is  not  well  known.  This  is  a  practical  problem, 
particularly  in  Arctic  environments  where  the  bathymetry  is  often  poorly  known  and  diSicult 
to  measure  due  to  ice  cover.  A  general  approach  to  the  problem  of  environmental  mismatch 
in  MFP  was  developed  by  Collins  and  Kuperman  [1]  based  on  generalizing  the  localization 
problem  to  include  uncertain  environmental  parameters  as  well  as  source  location  in  the 
search  for  an  optimal  match  with  the  acoustic  data.  Dosso  [2]  applied  this  method  to 
source  localization  in  environments  with  uncertain  bathymetry.  Here,  a  broadband  MFP 
algorithm  is  described  which  searches  simultaneously  for  the  location  (range  and  bearing) 
and  relative  strength  of  multiple  sources,  and  includes  corrections  to  the  bathymetry  along 
the  propagation  path  to  each  source  in  the  inversion.  The  approach  is  illustrated  with 
a  synthetic  example  representing  a  distribution  of  active  ice-ridge  building  events  in  the 
Arctic.  The  inversion  is  also  applied  to  a  set  of  Arctic  ambient  noise  measurements  collected 
in  the  Lincoln  Sea. 

1.  Inversion  Algorithm 

The  goal  of  MFP  is  to  determine,  the  location  of  one  or  more  sources  which  minimize 
the  mismatch  between  measured  and  modelled  acoustic  fields.  For  a  broadband  source,  the 
mismatch  E  can  be  defined  as 


E  =  \- 


>/u 


r"  P+,  \R,m 


u 


1/2 


>=1  ^QZ=:. 


L^p= 


E{“ 


q=p-j-l  ^f=fi 


-[1/2 


(1) 


where  H  is  the  number  of  hydrophones,  fi  and  /„  are  the  lower  and  upper  limits  of  the 
frequency  band,  and  Rpg  —  {RpR*g)  and  Mpg  =  Mp  M*  are  the  cross  spectra  of  the  measured 
and  modelled  acoustic  fields  Rp  and  Mp.  Equation  (1)  is  a  multiple-frequency  variant  of  the 
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Bartlett  processor  which  omits  the  autospectra  of  the  cross-spectral  matrices.  The  Bartlett 
processor  is  used  because  it  is  less  sensitive  to  environmental  mismatch  than  high-resolution 
processors.  Processing  multiple  frequencies  decreases  the  number  and  amplitude  of  sidelobes, 
since  these  tend  to  occur  at  different  ranges  for  different  frequencies.  The  modelled  acoustic 
field  is  determined  by  specifying  pressure  levels  at  the  receiver  as 


k{f)p{ri,Zi,9i,Zg,f) 

1 1/2  ’ 

Ef=i  \piri,Zu9i,ZjJ)\‘^] 

where  S  is  the  number  of  sources  and  (/)  is  the  total  power  received  at  the  array  for  source 
i.  The  acoustic  pressure  p  at  receiver  depth  Zg  due  to  a  source  at  is  calculated 

using  an  Nx2D  adiabatic  normal  mode  model 

exp(i  /o°  kj[d{r,  Oj 

/cTST)] 

where  N  is  the  number  of  modes  that  propagate  from  source  to  receiver,  (l>j[d{r,  9),  z]  is  the 
jth.  mode  function  at  depth  z  for  a  water  depth  d{r,  9),  and  kj[d{r,  0)]  is  the  jth  wavenumber. 
The  mode  functions  and  wavenumbers  were  pre-calculated  for  20-m  increments  in  water 
depth  and  stored  in  look-up  tables  for  fast  reference.  Within  the  invesion  algorithm,  linear 
interpolation  of  the  table  entries  was  then  used  to  determine  modal  values.  This  provides 
highly  efiicient  computation  of  the  matching  fields.  This  is  of  key  importance  to  the  multi¬ 
dimensional  inversion  developed  here,  which  requires  the  matching  fields  be  modelled  an 
extremely  large  number  of  times. 

The  matched-field  localization  algorithm  searches  simultaneously  for  the  range,  bearing 
and  relative  strength  of  a  specified  number  of  sources,  together  with  corrections  to  the 
bathymetry  model  along  the  propagation  path  to  each  source.  This  requires  an  extremely 
large  parameter  space  be  examined  to  find  the  parameter  set  which  minimizes  the  mismatch. 
The  parameter  search  is  carried  out  using  the  method  of  simulated  annealing  which  provides 
an  efficient  stochastic  search  algorithm  that  avoids  becoming  trapped  in  unfavourable  local 
minima  [1,  2]. 


Pin,Zi,9i,Zg)  =  (pj{d{ri,9i),Zi]  0j[d(0, 0), zj 

j=i 


)] 


dr 
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2.  Synthetic  Example 

Inversion  for  multiple  sources  with  uncertain  bathymetry  is  investigated  in  this  section 
using  a  synthetic  example  representing  a  distribution  of  active  ice-ridge  building  events  in  an 
Arctic  environment.  The  synthetic  case  is  designed  to  simulate  ambient  noise  measurements 
collected  in  the  Lincoln  Sea  which  are  analyzed  in  Sec.  3.  The  environment  and  array 
geometry  at  the  measurement  site  were  adopted  for  the  synthetic  study  and  are  summarized 
here.  The  hydrophone  array  was  located  in  430-m  of  water  and  consisted  of  vertical  and 
horizontal  sub-arrays  of  24  and  7  sensors,  respectively.  The  horizontal  sub-array  included  one 
off-axis  hydrophone  to  break  the  left-right  symmetry.  The  sound-speed  profile  was  measured 
at  the  array  site  and  a  number  of  experiments  were  carried  out  to  measure  the  ocean-bottom 
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Fig.  1  Mismatch  E  vs.  number  of  sources  S  for  the  synthetic  example. 


geoacoustic  properties  [3].  The  bathymetry  was  taken  from  the  chart  of  the  Arctic  basin 
published  by  the  Naval  Research  Laboratory  (NRL)  [4].  This  chart  was  digitized  to  400-km 
range  from  the  array  site  along  radial  paths  every  1°  in  bearing.  Synthetic  acoustic  data  and 
matching  fields  were  computed  for  the  frequency  band  10-30  Hz  in  5-Hz  increments  using 
the  adiabatic  normal  mode  model.  The  applicability  of  the  adiabatic  approximation  to 
this  environment  was  investigated  in  a  simulation  study  by  Zala  and  Ozard  [5],  who  found 
excellent  agreement  between  an  adiabatic  normal  mode  model  and  a  parabolic  equation 
model  at  20  Hz  for  slopes  <  3°.  The  largest  bottom  slopes  along  radial  paths  from  the 
array  site  are  <  2.4°;  hence,  the  adiabatic  normal  mode  model  should  be  applicable  to  this 
environment. 

The  number  of  sources  contributing  to  the  total  acoustic  field  may  be  estimated  by 
examining  the  mismatch  E  returned  by  the  inversion  algorithm  as  a  function  of  the  number 
of  sources  S  sought  in  the  inversion.  Figure  1  shows  E  sls  a  function  of  S  when  six  sources 
were  used  as  input  to  compute  synthetic  acoustic  data  (the  source  distribution  is  shown  in 
Fig.  2a).  The  search  interval  for  the  bathymetry  was  ±50  m.  In  Fig.  1,  the  mismatch  E 
decreases  as  the  number  of  sources  S  increases  to  six.  Increasing  S  beyond  six  does  not 
significantly  decrease  E.  Figure  2  shows  the  estimated  source  locations  and  relative  levels 
determined  by  the  inversion  algorithm.  When  the  number  of  sources  is  underestimated, 
the  stronger  sources  are  localized.  When  the  correct  number  of  sources  is  sought,  excellent 
estimates  of  the  source  locations  and  levels  are  obtained.  When  the  number  of  sources  is 
overestimated,  excellent  estimates  for  the  input  sources  are  still  obtained,  with  additional 
weak  sources  included  to  account  for  the  required  number.  These  additional  sources  do  not 
contribute  significantly  to  the  acoustic  field. 
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Fig.  2  True  source  distibution  (a,  shaded  circles)  and  distributions  determined  by  inversion 
for  2,  4,  6,  8  and  10  sources  in  (b)-(f).  Circle  size  indicates  relative  source  strength. 


3.  Arctic  Ambient  Noise  Inversions 

The  measured  data  analyzed  in  this  section  consist  of  six  two-minute  samples  of  the 
ambient  noise  recorded  in  the  Lincoln  Sea.  Five  of  the  samples,  denoted  samples  B-F,  were 
dominated  by  distant  ice-ridging  noise.  The  sixth,  sample  A,  was  dominated  by  local  thermal 
ice  cracking  and  was  included  in  the  analysis  for  comparative  purposes.  Matching  fields  for 
the  inversion  were  calculated  using  the  adiabatic  normal  mode  model  (Eq.  3).  An  important 
factor  infiuencing  propagation  in  the  Arctic  is  the  scattering  of  acoustic  energy  at  the  rough 
under-side  of  the  ice.  In  computing  the  matching  fields,  the  Burke-Twersky  ice-scattering 
model  [6]  was  employed,  which  represents  under-ice  ridges  by  a  random  distribution  of  ellip¬ 
tical  half-cylinders.  The  bathymetry  was  allowed  to  vary  in  the  inversions  by  ±50  m  from 
values  from  the  NRL  chart  [3]. 

Figure  3  shows  the  mismatch  £■  as  a  function  of  the  number  of  sources  S  sought  in  the 
inversion  for  the  six  ambient  noise  samples.  For  the  ice- ridging  samples  (B-F),  E  decreases 
rapidly  as  S  increases  to  about  6-8,  then  decreases  much  more  slowly.  For  seven  sources 
E  <  0.01  for  each  of  these  samples.  The  fit  to  the  data  is  illustrated  in  Fig.  4,  which 
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Fig.  3  Mismatch  E  vs.  number  of  sources  S  for  five  ambient  noise  sampes. 


shows  the  magnitude  of  the  cross-spectral  matrices  for  sample  E  and  for  the  modelled  source 
distribution  with  seven  sources.  All  the  major  features  of  the  measured  cross  spectra  are 
reproduced  by  the  modelled  source  distribution.  Thus,  6-8  sources  appear  to  be  sufficient  to 
model  the  ice- ridging  noise  samples  used  in  this  study.  For  sample  A,  the  mismatch  in  Fig.  3 
is  much  larger  than  for  the  ice-ridging  samples  and  does  not  decrease  significantly  as  the 
number  of  sources  is  increased  up  to  20.  This  indicates  that  the  inversion  method  does  not 
have  so  many  degrees  of  freedom  that  it  can  simply  reproduce  any  acoustic  measurements, 
such  as  local  (non-modal)  ice-cracking  noise. 

Each  of  the  inversions  represented  in  Fig.  3  involved  computing  approximately  7  x  10^ 
matching  fields,  and  required  ~2  hours  CPU  time  on  a  DEC  Alpha  computer  (computation 
speed:  20  Mflops).  To  verify  that  the  parameter-space  search  was  sufficient,  several  inversions 
were  repeated  with  the  number  of  matching  fields  increased  by  a  factor  of  ten,  with  no 
appreciable  change  in  the  mismatch.  To  verify  that  the  ±50  m  search  interval  for  the 
bathymetry  was  sufficient,  several  cases  were  repeated  with  a  search  interval  of  ±250  m, 
with  no  change  in  mismatch. 

The  source  distributions  determined  by  the  inversions  for  seven  sources  are  shown  in 
Fig.  5.  This  figure  also  includes  the  results  of  broadband  maximum-likelihood  plane- wave 
beamforming  applied  to  the  acoustic  data  recorded  on  the  horizontal  array  for  the  frequency 
band  10-30  Hz.  The  results  suggest  that  sources  located  in  the  deep  water  to  the  north  of 
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Fig.  4  Magnitude  of  cross-spectral  matrices  for  measured  (a,  c,  e)  and  modelled  (b,  d,f) 
fields.  Frequencies  are  10  Hz  (a,  b),  20  Hz  (c,  e)  and  30  Hz  (e,f). 

the  recording  site  dominate  the  ice-ridging  samples  B-F.  The  ice-cracking  sample,  sample  A, 
has  much  flatter  directivity  than  the  ice-ridging  samples,  indicating  an  essentially  isotropic 
distribution  of  sources. 

Given  that  inverse  problems  are  inherently  non-unique,  the  good  agreement  between 
measured  and  modelled  fields  in  Fig.  4  does  not  necessarily  indicate  that  the  sources  are 
correctly  localized.  A  test  of  the  result  is  to  obtain  consistent  source  distributions  for  different 
initializations  of  the  inversion  algorithm.  Unfortunately,  this  has  not  proved  to  be  the  case  for 
the  measured  ambient  noise  data.  Source  ranges  were  found  to  be  ambiguous,  i.e.,  different 
initializations  of  the  inversion  lead  to  source  distributions  with  similar  mismatches  that 
differed  significantly  in  source  ranges.  However,  source  bearing  estimates  were  found  to  be 
consistent  with  respect  to  different  initialization  of  the  inversion.  In  addition,  re-initializing 
the  inversion  algorithm  for  different  numbers  of  sources  did  not  change  the  conclusion  that 
6-8  sources  were  required  to  model  the  ice-ridging  noise. 

There  are  a  number  of  possible  reasons  for  the  ambiguity  in  range  for  the  measured  data. 
These  include  inadequacies  of  the  ice-scattering  model,  mode  coupling  and  three-dimensional 
propagation  effects,  errors  in  the  location  of  the  hydrophones,  non-modal  noise  from  other 
sources,  and  uncertainty  in  environmental  parameters  other  than  bathymetry. 
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Fig.  5  Source  distributions  (circles)  for  samples  A-F  in  (a)-(f),  respectively.  Dashed  lines 
indicate  beamforming  results. 

4.  Summary 

A  matched-field  inversion  algorithm  is  described  for  localizing  multiple  sources  with  a 
limited  knowledge  of  bathymetry.  The  method  provides  excellent  results  for  synthetic  ex¬ 
amples  representing  a  distribution  of  ice-ridging  noise  sources  in  the  Arctic.  Application  to 
ambient  noise  measurements  provided  source  distributions  which  reproduce  the  measured 
fields  and  are  unambiguous  in  terms  of  the  number  and  bearings  of  the  sources;  however, 
source  ranges  could  not  be  estimated  unambiguously. 
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Abstract 

Matched  field  processing  of  broadband  SUS  charge  data  was  performed  to  evaluate  a  coherent  broad¬ 
band  object  function.  This  object  function  includes  absolute  phase  and  amplitude  information  at  all 
frequencies,  thereby  matching  the  shapes  and  absolute  timing  of  pressure  waveforms  measured  on  a  ver¬ 
tical  array.  The  coherent  broadband  object  function  is  implemented  together  with  the  classical  Bartlett 
object  function,  and  an  incoherent  broadband  function  obtained  by  averaging  the  Bartlett  function  over 
frequency.  All  three  functions  are  applied  to  determine  their  effectivenesses  for  inverting  experimental 
geometry  and  ocean  bottom  parameters. 


Introduction 

Matched  Field  Processing  (MFP)  has  been  discussed  in  recent  underwater  acoustics  literature.  The  method 
is  applied  as  follows:  certain  characteristics  of  the  pressure  field  measured  on  an  array  of  hydrophones  are 
matched  with  those  predicted  by  a  numerical  model  by  varying  the  inputs  to  the  model.  The  inputs  causing  the 
best  match  are  taken  as  estimates  of  the  corresponding  parameters  present  in  the  vicinity  of  the  experiment. 

The  inputs  to  the  numerical  model  normally  include  the  positions  of  the  source  and  the  receiver(s),  and 
a  set  of  parameters  describing  the  geoacoustic  environment.  We  refer  to  a  specific  set  of  inputs  to  the 
propagation  model  as  an  environmental  model  m.  It  is  denoted  in  bold  to  indicate  that  it  represents  the 
values  of  several  parameters  at  once.  An  object  function  E{m.)  is  defined  to  indicate  the  degree  of  match 
for  specific  models.  The  problem  then  is  to  maximize  the  object  function  over  a  search  space  defined  by  the 
possible  ranges  for  each  of  the  components  of  m. 

The  object  or  cost  function  dictates  which  information  in  the  pressure  field  is  matched.  The  following 
sections  describe  possible  forms  of  the  object  function.  The  first  type  discussed  is  the  single  frequency 
Bartlett  object  function.  Next,  extension  of  the  single  frequency  function  to  the  incoherent  broadband  object 
function  is  described.  Finally  the  coherent  broadband  object  function  is  developed,  and  discussed  in  relation 
to  the  incoherent  cases. 


I  Experiment 

Experimental  data  was  aquired  by  the  Defence  Research  Establishment  Pacific  (DREP)  during  the  Pacific 
Shelf  sea  trials  in  1993  off  the  coast  of  Vancouver  Island.  The  shot  data  processed  for  this  work  were  recorded 
during  a  shot  run  in  which  a  series  of  0.82  kg  SUS  charges  were  dropped  at  increasing  range  from  a  vertical 
array  which  was  stationary  at  approximately  49.02®N  126.86°W,  where  the  water  depth  was  approximately 
400  m.  The  charges,  set  to  detonate  at  200  m,  were  dropped  at  ranges  between  0.7  km  and  12  km  from 
the  array,  along  a  course  which  approximately  paralleled  the  400  m  depth  contour.  The  array  consisted  of 
16  hydrophones  equally  spaced  with  15  m  separations,  with  the  depth  of  the  top  hydrophone  at  80  ±  5  m 
throughout  the  run.  Hydrophone  pressure  was  sampled  continuously  at  a  rate  of  1500  samples  per  second 
on  all  hydrophones.  These  data  were  transmitted  via  a  FM  link  from  the  array  to  a  monitoring  ship,  where 
they  were  written  directly  to  exabyte  tape  files. 


II  Object  Functions 


II. 1  Incoherent  object  functions 


The  most  common  formulation  of  the  object  function  in  MFP  is  the  single  frequency  classical  Bartlett  power 
function: 


E<3i(m)  Pi 
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where  N  is  the  number  of  hydrophones,  and  Q2(m)  and  Pi  are  the  predicted  and  measured  complex  Fourier 
components  at  one  frequency  on  the  z'th  hydrophone,  and  *  represents  complex  conjugation.  jE(m)  ranges 
from  0  to  +1,  with  -fl  representing  a  perfect  match.  This  occurs  when  Q2(m)  =  a  Pi  for  all  i,  where  a  is 
any  complex  constant. 

An  incoherent  broadband  object  function  can  be  obtained  from  a  weighted  average  of  equation  1  over 
frequency.  If  M  equally  weighted  frequencies  are  included  the  multiple  frequency  object  function  is 
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II. 2  Coherent  Object  Functions  (Waveform  matching) 

A  method  for  full  field  coherent  broadband  inversion  was  suggested  by  researchers  in  the  geophysics  field  [1]. 
This  method  involves  matching  synthetic  waveforms  qi{t),  referred  to  as  synthetic  seismograms,  with  measured 
waveforms  pi{i).  Prediction  of  the  synthetic  waveforms  at  each  of  N  hydrophones  requires  knowledge  of  the 
source  waveform  s(f),  and  the  waveguide's  impulse  response  functions  gi{t),  i  =  Synthetic  waveforms 

are  obtained  by  convolving  the  source  waveform  by  the  impulse  response  functions: 

qi{t)  =  s{t)  *  gi{t).  (3) 

The  convolution  can  usually  be  done  more  efficiently  in  the  frequency  domain.  By  the  convolution  theorem, 
the  Fourier  transform  of  equation  3  can  be  written. 

Qi{u>)  =  S{u)Gi{u!),  (4) 
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where  Qi{u))^  S{(jj)^  and  Gt(a;)  are  the  Fourier  transforms  of  qi[t)^  s{t)j  and  gi{t)  respectively.  A  numerical 
model  is  used  to  predict  Gi{uj)  for  each  environmental  model  m.  The  procedure  for  calculating  qi{t)  then  is 
to  take  the  inverse  Fourier  transform  of  the  right  side  of  equation  4:  qi{t)  =  (5(w)Gj(a;)). 

Because  we  wish  to  match  waveforms,  a  good  choice  for  the  object  fxmction  is  the  correlation  between 
synthetic  and  measured  waveforms.  In  order  to  be  consistent  with  the  previously  described  object  functions, 
we  define  a  new  object  function  as  the  normalized  sum  over  the  hydrophones  of  the  squared  correlations  of 
measured  and  predicted  time  pressure  series: 
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Application  of  equation  5  requires  that  the  absolute  timing  of  the  measured  waveforms,  referred  to  the  time  of 
source  radiation,  be  known.  Absolute  timing  is  unknown  unless  the  source  signal  is  monitored  at  the  source 
location,  and  the  absolute  time  reference  between  this  and  the  signals  recorded  at  the  array  is  established. 
A  method  for  treating  this  problem,  when  absolute  timing  is  not  available,  is  to  search  for  the  time  delay 
r  which  maximizes  the  correlation  between  qi{t  +  r)  and  Pi{t).  This  is  done  by  replacing  the  correlation 
coelficient  in  the  numerator  of  equation  5  by  the  correlation  function.  The  object  function  is  then  calculated 
from  the  maximum  of  the  correlation  function: 
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The  search  for  optimum  time  delay  is  performed  outside  of  the  summation  over  hydrophones  because  the 
same  delay  is  common  to  all  hydrophones.  This  object  fimction  can  be  calculated  in  the  frequency  domain 
making  use  of  the  correlation  and  ParsevaFs  theorems: 
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The  integrals  over  continuous  frequency  can  be  replaced  by  summation  over  discrete  frequencies  if  the  source 
signal  is  band  limited.  Similarly,  if  the  signal  contains  discreet  tonals,  the  summations  can  be  limited  to  those 
frequencies  only.  In  either  case,  with  M  discrete  frequencies,  equation  7  becomes 
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An  interesting  special  case  of  equation  8  is  when  there  is  only  a  single  frequency  ujq.  For  this  case,  equation 
8  reduces  to 
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which  is  exactly  equivalent  to  the  Bartlett  object  function  of  equation  1.  The  absolute  time  uncertainty  which 
required  searching  over  time  delays  does  not  atfect  this  last  expression. 


Ill  Application  of  MFP 

III.l  Environment  and  geometry  models 

The  sediment  at  the  seafloor  was  fine  to  medium  grain  sand.  A  standard  analysis  of  strong  head  wave  arrivals 
present  in  the  data  indicated  the  existence  of  at  least  two  sub-bottom  layers  having  compressional  speeds  of 
1890±30  m/s  and  2070±30  m/s.  The  top  of  the  first  of  these  layers  is  at  least  150  m  beneath  the  seafloor, 
and  the  thickness  of  this  layer  is  approximately  200  m.  The  geoacoustic  environment  chosen  for  inversion 
modeling  purposes  consisted  of  three  homogeneous  elastic  ocean  bottom  layers.  Several  of  the  parameters 
describing  these  layers  were  fixed  at  values  obtained  by  matching  the  above  layer  compressional  speeds  with 
those  from  similar  sediments  categorized  in  Hamilton’s  compilation  of  seafloor  sediment  properties  [2].  In 
reference  to  Hamilton’s  work,  the  sediment  types  attributed  to  the  three  layers  were  sand,  gravel  and  morain. 
The  fixed  parameter  values  assigned  to  these  layers  are  given  in  table  1. 

Models  representing  the  geoacoustic  and  geometry  information  to  be  inverted  for  were  described  by: 
m  =  {  dsj  T),  ^1,  t2,  Cl,  C2j  C3},  where  dg  is  the  source  depth,  dr  is  the  depth  of  the  top  hydrophone,  r 

is  the  horizontal  source  to  array  range,  D  is  the  bottom  depth,  ti  and  t2  are  the  thicknesses  of  the  top  two 
layers,  and  ci,  C2,  and  C3  are  the  compressional  speeds  in  the  layers.  Additionally,  the  true  measured  ocean 
sound  speed  profile  was  used,  and  the  measured  array  tilt  angle  was  fixed  at  5®. 

III.  2  Inversion 

Because  computational  time  is  the  major  limiting  factor  for  broadband  propagation  calculations,  the  normal 
mode  model  ORCA  [3]  was  chosen  as  the  forward  model.  ORCA  eflficiently  calculates  broadband  waveguide 
transfer  fxmctions  by  tracking  mode  eigenvalues  through  a  range  of  frequencies.  Increasing  the  specified 
frequency  resolution  of  the  transfer  function  causes  very  little  change  in  computation  time.  This  characteristic 
made  it  possible  and  eflficient  to  apply  the  object  functions  at  a  large  number  of  frequencies.  For  our  broadband 
trials,  136  equally  spaced  frequencies  between  10  Hz  and  60  Hz  were  used.  The  inversion  was  performed  on 
data  from  a  single  shot  at  approximate  range  3.2  km.  Execution  time  was  approximately  90  s  per  model  on 
a  150  MHz  R4000  Silicon  Graphics  workstation.  However,  this  time  included  calculation  of  all  ranges  (r)  at 
once. 

The  inversion  was  performed  using  a  grid  search  in  which  the  incoherent  broadband  object  function 
(equation  2)  and  the  coherent  broadband  object  function  (equation  8)  were  evaluated  for  a  set  models  defined 
by  discrete  values  within  specified  limits  for  each  of  the  components  of  m.  Additionally,  the  single  frequency 
object  function  (equation  1)  was  evaluated  at  a  frequency  of  50  Hz  for  this  same  set  of  models.  The  discrete 
parameter  values  for  m  are  defined  by  the  “Minimum” ,  “Maximum” ,  and  “Num  Steps”  entries  in  table  2.  On 


Layer 

Shear  speed 

Density 

P-atten 

S-atten 

1.5  g/cm^ 

0.8  dB/A 

2.0  dB/A 

1.6  g/cm^ 

1.5  dB/A 

1.7  g/cm^ 

1.0  dB/A 

Table  1:  Fixed  parameter  values  for  bottom  layers. 
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Parameter 

Minimum 

Maximum 

Num  Steps 

Inversion  Result 

r 

3.1  km 

3.3  km 

5 

3.25  km 

D 

425  m 

3 

425  m 

dj' 

85  m 

3 

75  m 

240  m 

3 

h 

240  m 

3 

Cl 

1700  m/s 

7 

C2 

2000  m/s 

4 

C3 

fSsSaBBam 

2300  m/s 

5 

Table  2:  Parameter  search  ranges  and  coherent  broadband  results. 


each  iteration,  ORCA  generates  the  broadband  output  transfer  function  Gf(m),  A  predicted  field  spectrum 
Qj(m)  is  calculated  according  to  equation  4.  The  source  spectrum  5(a;),  used  for  this  calculation,  was 
obtained  from  high  resolution  SUS  charge  recordings  collected  during  a  controlled  experiment  [4]. 

The  model  which  produced  the  largest  broadband  coherent  object  function  value  was  determined.  The 
parameter  values  for  this  model  are  given  in  table  2.  As  mentioned  in  section  11.2,  application  of  the  broadband 
coherent  object  function  is  equivalent  to  waveform  matching.  To  demonstrate  this,  the  measured  waveforms 
and  predicted  waveforms  based  on  the  best  model  are  shown  together  in  figure  1 . 


Figure  1:  Modeled  and  predicted  waveforms 


III.3  Sensitivity  Study 

In  order  to  obtain  a  measure  of  the  robustness  of  the  individual  parameter  estimates,  we  observed  the  amounts 
of  change  in  the  object  functions’  values  as  the  parameters  were  varied  through  their  respective  search  ranges. 
This  was  performed  for  all  parameters  of  m  with  the  exception  of  source  depth,  which  was  fixed  at  196  m. 
When  each  parameter  was  varied,  the  remaining  parameters  were  fixed  at  the  optimal  values  given  in  table  2. 
The  results  of  this  study  indicated  that  varying  the  parameters:  source  to  receiver  range,  bottom  depth, 
and  array  depth,  causes  significant  change  in  all  object  functions.  The  broadband  functions  are  particularly 
sensitive  to  the  upper  layer  compressional  speed,  in  contrast  to  the  single  frequency  object  function  which 
is  fairly  insensitive  to  this  parameter,  and  shows  two  local  maxima  over  its  search  range.  All  three  object 
functions  are  very  insensitive  to  the  remainder  of  the  ocean  bottom  parameters.  An  interesting  feature  is  that 
the  single  frequency  and  broadband  incoherent  functions  have  consistently  and  significantly  higher  values 
than  the  broadband  coherent  function.  The  reason  for  this  feature  is  unknown. 
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III. 4  Source  Localization  Problem 


The  sensitivity  study  indicated  that  all  three  object  functions  were  relatively  insensitive  to  the  geoacoustic 
parameters  of  the  lower  two  subbottom  layers,  but  that  they  were  sensitive  to  experimental  geometry.  Con¬ 
sequently,  the  methods  are  expected  to  perform  well  for  the  common  problem  of  localizing  a  source.  A  test 
of  MFP  source  localization  was  performed  by  fixing  all  geo  acoustic  parameters  at  the  optimal  values  shown 
in  table  2.  The  object  functions  were  evaluated  as  a  function  of  source  depth  and  range  from  the  array.  The 
localizations  tested  a  grid  of  81  ranges,  between  2.0  km  and  4.0  km,  by  41  depths,  between  10  m,  and  410  m. 
Contour  plots  of  object  function  versus  source  range  and  depth  are  referred  to  as  ambiguity  surfaces.  Figures 
2,  3,  and  4  show  the  ambiguity  surfaces  respectively  for  the  single  frequency,  incoherent  broadband  and 
coherent  broadband  object  functions.  All  three  processors  give  very  good  results,  with  strong  peaks  at  the 
source  position:  3.25  km  range  and  196  m  depth.  The  single  frequency  object  function  performs  suprisingly 
well  even  though  it  uses  significantly  less  acoustic  field  information.  In  fact,  the  peak  at  the  source  location  is 
even  sharper  than  that  of  both  broadband  functions.  This  is  likely  due  to  the  fact  that  the  broadband  object 
functions  included  lower  frequencies  for  which  there  is  significant  energy  in  the  source  spectrum.  These  lower 
frequencies  have  broader  main  lobes  in  the  ambiguity  surface. 

Secondary  lobes  for  the  single  frequency  processor  are  distributed  fairly  evenly  over  the  ambiguity  surface. 
Both  the  broadband  coherent  and  broadband  incoherent  processors  have  fewer,  and  smaller  secondary  lobes 
than  the  single  frequency  processor.  Additionally,  these  lobes  are  elongated,  showing  some  coupling  between 
the  range  and  depth  parameters. 


Figure  2:  Single  frequency  Bartlett  ambiguity  surface 


Figure  3:  Incoherent  Broadband  ambiguity  surface 
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Figure  4:  Coherent  Broadband  ambiguity  surface 


Conclusions 

Matched  Field  Processing  of  transient  SUS  charge  data  was  used  to  evaluate  and  compare  coherent  broad- 
band,  incoherent  broadband,  and  single  frequency  object  functions  for  the  problem  of  geoacoustic  parameter 
inversion.  The  evaluations  were  performed  on  data  from  a  charge  at  3.2  km,  in  a  frequency  band  between 
10  Hz  and  60  Hz  for  the  broadband  processors,  and  at  50  Hz  for  the  single  frequency  processor.  The  results  of 
these  evaluations  indicate  that  all  three  processors  are  sensitive  to  experimental  geometry,  but  none  was  found 
to  be  sensitive  to  the  geoacoustic  parameters  for  the  second  and  third  bottom  layers.  Only  the  broadband 
processors  are  sensitive  to  the  compressional  speed  in  the  top  ocean  bottom  layer. 

An  important  feature  of  the  coherent  broadband  processor,  applied  in  this  evaluation,  is  that  it  required 
searching  over  time  delays.  This  step  was  necessary  because  the  absolute  timing  of  the  measured  signals  with 
respect  to  source  detonation  was  unavailable.  If  the  source  had  been  monitored  and  recorded  with  an  absolute 
time  reference  to  the  signals  on  the  array,  searching  over  time  delay  would  not  be  necessary  and  calculation 
of  the  full  correlation  function  would  not  be  required.  More  importantly,  absolute  propagation  times  as  well 
as  hydrophone  to  hydrophone  differential  times  for  all  acoustic  paths  would  be  matched.  Greater  sensitivity 
to  parameters  affecting  absolute  propagation  times,  such  as:  range,  ocean  depth,  and  bottom  layer  thickness 
to  layer  speed  ratios,  would  be  expected. 
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Abstract 

The  determination  of  shallow  water  bottom  properties  can  be  extremely  difficult,  but  among  the  most 
promising  methods  are  remote  sensing  techniques  which  exploit  the  sensitivity  of  propagating  acoustic 
fields  to  such  properties.  In  particular,  there  are  the  Matched  Field  Processing  (MFP)  based  inverse 
techniques  which  consider  both  the  acoustic  phases  and  amplitudes  in  their  hunt  for  estimates  of  such 
environmental  properties  as  sediment  layer  thicknesses,  densities,  and  sound-speeds.  This  paper  focusses 
primarily  on  recent  efforts  by  the  author  to  estimate  such  parameters  for  two  experimental  data  sets 
(Mediterranean  Sea  north  of  Elba  and  Hudson  Canyon).  However,  an  important  but  unresolved  issue 
for  all  the  inverse  techniques  (not  just  the  MFP  based  ones)  remains  that  of  uniquely  determining  the 
“true”  or  “solution”  set  of  parameter  values. 


Introduction 

Matched  Field  Processing  (MFP)  has  recently  been  the  explicit  basis  for  a  variety  of  methods  for  the 
estimation  of  shallow  water  bottom  properties.  The  particular  MFP  approach  to  be  discussed  here  is  similar 
to  the  MFP  tomographic  method  used  for  the  estimation  of  deep  ocean  sound-speed  profiles,  except  that  now 
more  parameters  need  to  be  determined  (four  to  six  per  sediment  layer  excluding  elasticity  considerations) , 
the  frequencies  to  be  processed  need  to  be  higher  (50  to  500Hz) ,  the  source  to  array  ranges  are  much  shorter 
(less  than  20km) ,  and  the  field  sensitivities  to  the  individual  parameters  vary  by  scenario. 

In  this  paper  we  process  and  analyze  experimental  test  data  from  two  shallow  water  scenarios:  one  in  the 
Mediterranean  Sea  north  of  Elba  and  one  in  the  Hudson  Canyon  off  the  New  Jersey  Continental  shelf.  For 
each  we  assume  no  range  or  azimuthal  variability,  a  single  source,  and  a  single  nearly  vertical  array.  We  use 
the  KRAKEN  normal  mode  model[2]  and  concentrate  on  the  performance  of  the  Minimum  Variance  (MV) 
processor[4]. 

I  Discussion 

The  Mediterranean  Sea  data  involved  shot  sources  which  produced  a  broadband  of  frequencies  (from  60- 
420Hz)  displaying  strong  components  for  frequencies  around  lOOHz  over  the  entire  array  (we  selected  a 
uniform  subarray  of  9  phones  at  8m  spacings  with  the  top  phone  nominally  at  35m  depth).  The  processing 
and  inversion  were  performed  at  100.6Hz  and  eventually  assumed  two  sediment  layers  over  a  half-space  (for 
a  total  of  eight  environmental  parameters  to  be  determined).  This  resulted  in  several  very  good  MV  peak 
values  (approximately  0.5  out  of  a  maximum  possible  of  1.0)  near  the  “known”  source  location  (additionally 
refined  by  the  processing).  The  MV  based  method  was  seen  to  be  quite  sensitive  to  phone  location  errors, 
and  so,  an  iterative  algorithm  was  introduced  to  allow  for  the  inclusion  of  phone  coordinates  as  parameters. 
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Subsequently,  phone  range  displacements  of  less  than  1,0m  resulted  in  improved  MV  peak  values  up  to  0.7. 
Complete  details  can  be  found  in  Tolstoy  (1996).  Unfortunately,  these  new  parameters  do  not  result  in 
such  strong  MFP  values  at  neighboring  frequencies,  suggesting  that  the  “true”  set  of  values  has  yet  to  be 
determined.  It  should  be  noted  that  the  Linear  processor  peaks  were  frequently  near  the  maximum  value 
1.0.  The  to-be-expected  tradeoff  is  that  while  an  MV  based  method  is  far  more  sensitive  to  measurement 
errors,  it  can  more  easilty  suppress  sidelobe  solutions.  It  seems  that  the  Linear  processor  has  such  an  easy 
time  finding  peaks  near  1.0  that  only  a  high  resolution  method  may  have  sufficient  sidelobe  suppression  to 
find  a  “true”  solution.  Perhaps  we  can  generally  conclude  that  if  an  inversion  method  produces  estimates 
which  do  not  result  in  strong  MV  peak  values,  then  we  have  a  good  indication  of  the  presence  of  some  kind 
of  unknown  mismatch,  and  the  inversion  results  should  be  suspect. 

Some  of  the  Hudson  Canyon  data  have  also  been  processed  through  the  MFP  MV  inversion  software 
where  these  data  were  selected  because  they  are  extremely  weU  documented[l],  because  they  are  readily 
available  by  anonymous  ftp,  and  because  ear  her  inversion  efforts  have  already  made  excellent  progress  in 
the  broadband  processing  of  those  data  for  the  purposes  of  source  tracking[3].  For  our  processing  we  used 
a  uniform  subarray  (an  8  phone  subset  at  7.5m  spacings  with  the  top  phone  at  19.95m)  and  selected  the 
50Hz  CW  signal  at  a  nominal  range  of  4.58km.  After  the  inversion,  we  arrived  at  a  six  parameter  set  (single 
sediment  layer  over  a  half  space)  which  produced  an  MV  peak  value  of  0.76  (and  a  number  of  sets  producing 
values  near  0.7  with  Linear  processor  peaks  at  1.00).  After  including  iterations  for  phone  displacement  we 
finally  arrived  at  a  top  MV  peak  value  of  0.83.  More  details  can  be  found  in  Tolstoy  and  Michalopoulou 
(1996).  Additional  frequencies  remain  to  be  examined.  However,  it  is  not  hkely  that  we  have  yet  found  the 
“true”  solution  set  for  these  data. 

At  this  point  it  seemed  like  a  good  idea  to  pause  rather  than  blindly  continue  the  searches  for  more 
parameter  sets.  The  MV  based  MFP  method  discussed  here  requires  progressive  searching  through  parameter 
spaces  which  are  enormous  (4  to  6  parameters  per  layer  each  with  dozens  to  hundreds  of  possible  discrete 
values,  plus  simultaneous  refinements  for  source  range  and  depth,  and  iterative  refinements  for  array  phone 
coordinates).  Moreover,  the  spaces  axe  very  non-convex  making  the  search  for  a  global  maximum  exceedingly 
difficult.  Unfortunately,  even  for  these  very  large  search  spaces,  one  major  uncertainty  (for  both  these  data 
sets)  remains:  have  we  modeled  each  situation  appropriately,  i.e.,  does  each  search  space  encompass  the 
“true”  environment?  Or  do  we  need  to: 

•  search  over  larger  intervals  for  our  parameter  values? 

•  introduce  more  environmental  variables,  e.g.,  more  sediment  layers? 

•  consider  more  errors,  e.g.,  in  the  assumed  water  column  sound-speed  profile? 

•  consider  range  dependence? 


Conclusions 

When  applied  to  the  two  sets  of  experimental  test  data  discussed  above,  an  important  result  of  the  single 
frequency  inversion  is: 

•  strong  processing  power  peaks  for  the  highly  sensitive  MV  processor,  i.e.,  peaks  greater  than  0.5. 

The  additional  considerations  of  refined  source  location  plus  array  deformation  (non-uniform  range  shifts) 
result  in  even  larger  MV  processing  peaks.  The  corresponding  Linear  processor  showed  MFP  peak  values  of 
1.00.  However,  the  strong  MV  peaks  did  not  generalize  across  neighboring  frequencies  suggesting  that  the 
“true”  parameters  have  not  yet  been  found. 

A  number  of  critical  issues  remain  unresolved.  In  particular: 

•  How  realistic  are  the  environmental  parameters  estimated  by  this  (or  any)  MFP  technique?  How  close 
are  they  to  the  “true”  values? 
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•  No  measurement  technique  offers  sufficient  resolution  of  the  bottom  parameter  values  with  respect  to 
the  MV  MFP  method  discussed  here.  So,  how  do  we  validate  the  inversion  results?  Will  it  be  necessary 
and/or  sufficient  to  achieve  results  with  high  MV  peak  values  (greater  than  0.5  out  of  a  maxmimum 
of  1.0)  over  a  range  of  frequencies? 

•  How  valid  are  the  suggested  array  deformations?  Simulations  under  ideal  conditions  converge  to  ‘‘true” 
array  shapes,  but  this  is  not  the  case  when  errors  (such  as  those  due  to  environmental  mismatch)  are 
present. 

•  When  do  we  stop  the  search? 
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Abstract 

Accurate  knowledge  of  array  shape  is  essential  for  carrying  out  full  wavefield  (matched  field)  processing. 
Direct  approaches  to  array  element  localization  (AEL)  include  both  non-acoustic  (tilt-heading  sensors)  and 
acoustic  (high  frequency,  transponder-based  navigation)  methods.  The  low  frequency  signature  emitted  from  a 
distant  source  also  can  be  used  in  an  inversion  approach  to  determine  array  shape.  The  focus  of  this  paper  is  on  a 
comparison  of  the  array  shape  results  from  these  three  different  methods  using  data  from  a  120  m  aperture  vertical 
array  deployed  during  SWellEx-3  (Shallow  Water  evaluation  cell  Experiment  #3).  Located  2  m  above  the 
shallowest  array  element  was  a  self-recording  package  equipped  with  depth,  tilt,  and  direction-of-tilt  sensors, 
thereby  permitting  AEL  to  be  performed  non-acoustically.  Direct  AEL  also  was  performed  acoustically  by 
making  use  of  transponder  pings  (in  the  vicinity  of  12  kHz)  received  by  high  frequency  hydrophones  spaced  every 
7.5  m  along  the  vertical  array.  In  addition  to  these  direct  approaches,  AEL  was  carred  out  using  an  inversion 
technique  where  matched-field  processing  was  performed  on  a  multi-tone  (50-200  Hz),  acoustic  source  at  various 
ranges  and  azimuths  from  the  array.  As  shown,  the  time-evolving  array  shape  estimates  generated  by  all  three 
AEL  methods  provide  a  consistent  picture  of  array  motion  throughout  the  6  hour  period  analyzed. 

[Work  supported  by  ONR,  Code  321  US.] 
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Abstract 

This  simulation  study  examined  several  optimization-based  approaches  to  inversion  for  geoacoustic 
parameters.  Using  the  program  Orca  to  compute  normal  mode  modek  for  realistic  range-independent 
waveguides,  the  sensitivities  of  and  the  dependencies  between  several  of  the  geoacoustic  parameters  were 
examined.  Two-dimensional  plots  of  the  normalized  Bartlett  processor  as  a  function  of  parameter  pairs 
were  often  multimodal  and  of  complex  structure,  and  certain  parameter  pairs  exhibited  interdependence. 

Using  a  rapidly  evaluated  model  function  with  similar  overall  characteristics  (i.e.,  multiple  local  optima, 
different  parameter  sensitivities,  and  parameter  interdependence),  a  comparison  of  inversion  algorithms 
was  done;  these  included  simulated  annealing,  a  genetic  algorithm,  and  a  combined  random  search/local 
optimization  approach.  The  performance  of  several  local  optimization  algorithms  was  also  compared 
using  this  model  function.  Based  on  the  residts,  a  strategy  was  adopted  which  consisted  of  an  initial 
random  search  stage  followed  by  local  optimization  (using  a  quasi- Newton  optimizer)  of  those  best 
matches  found  during  the  search.  This  two-stage  technique  was  implemented  in  an  MFP/MFI  software 
system  written  in  IDL. 


I  Introduction 

Successful  MFP  requires  a  knowledge  of  the  bathymetry  and  the  existence  of  a  suitable  geoacoustic  model 
for  the  sound  speed  profile  and  the  bottom  characteristics.  In  practice,  these  parameters  are  known  only 
approximately,  and  it  is  desirable  to  obtain  more  accurate  estimates  to  improve  the  effectiveness  of  MFP. 
This  aim  may  be  achieved  through  matched-field  inversion  (MFI)  of  data  obtained  using  a  source  at  a  se¬ 
ries  of  known  locations  in  the  environment.  Development  of  effective  inversion  procedures  requires  the  use 
of  appropriate  numerical  optimization  techniques  suited  to  the  nature  of  the  multimodal  multidimensional 
matching  functional,  or  ambiguity  surface,  to  be  optimized. 


II  Pairwise  Plots 

Using  a  vertical  20-element  array,  and  the  normal  mode  program  Orca  [1]  to  model  the  acoustic  fields,  the 
relations  between  various  parameter  pairs  of  geoacoustic  models  typical  of  the  Pacific  and  Arctic  continental 
shelves  were  examined.  An  example  of  these  results,  showing  the  normalized  Bartlett  processor  output  as 
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a  function  of  channel  depth  and  compressional  sound  speed  in  the  sediment  layer,  is  shown  in  Figure  1. 
This  and  other  plots  indicated  the  multimodal  nature  of  the  function  to  be  optimized,  the  wide  range 
of  sensitivities  of  the  various  parameters,  and  the  existence  of  complex  interdependencies  between  certain 
parameters. 
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Figure  1:  2-D  plot  of  normalized  Bartlett  processor  output  as  a  function  of  channel  depth  and  compressional 
sound  speed.  The  true  values  (100  m  and  1700  m/s)  are  in  the  centre  of  the  plot. 

Ill  Comparison  of  Optimization  Techniques  using  Model 
Functions 

To  fully  compare  the  performance  of  various  algorithms  for  inversion  using  a  geoacoustic  modelling  program 
would  be  extremely  time-consuming.  Instead,  a  model  function  with  the  properties  described  above  was 
designed  for  use  in  this  comparison.  The  form  of  this  function  was 

^  S  cos2[27f{(L(x  -  Xq)  ■  ffc)}], 

^  ib=l 

where 

X  is  an  n-element  vector  of  parameters, 
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xo  is  a  vector  of  true  values  for  the  n  parameters, 

K  is  the  number  of  cosine  terms, 

fjfc  is  the  Arth  vector  of  frequencies,  with  one  element  of  fk  per  parameter,  and 

L  is  the  lower  triangle  of  the  Cholesky  decomposition  of  a  ’parameter  interdependency  matrix  R,  with 
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This  function  has  a  global  maximum  of  1.0  when  x  =  xq-  Parameter  sensitivity  was  simulated  in 
each  run  by  taking  the  frequency  values  fkj  corresponding  to  the  jih  parameter  as  realizations  from  a 
uniform  distribution  on  [-Uj^-^Uj].  Increasing  uj  results  in  an  increase  in  sensitivity.  Typically,  four  to 
six  parameters  were  used,  with  a  factor  of  two  or  three  between  the  successive  values  of  Uj,  to  give  a  wide 
range  of  sensitivities.  The  bounds  for  the  parameters  were  set  to  give  approximately  10-20  peaks  in  the  2-0 
surfaces  computed  for  the  two  most  sensitive  parameters. 

Using  a  large  number  of  realizations  of  this  functional  form  (i.e.,  a  different  set  of  values  of  the  for  each 
run),  simulated  annealing  [2]  and  a  genetic  algorithm  [3]  were  compared  for  their  rate  of  reduction  in  the 
function  value,  convergence  of  the  parameters,  and  success  in  localizing  the  parameters  to  the  region  of  the 
true  global  optimum.  It  was  found  that  the  performance  of  simulated  annealing  and  the  genetic  algorithm 
was  quite  similar,  but  their  rates  of  convergence  were  very  slow.  Furthermore,  their  success  in  identifying  a 
point  in  the  region  of  the  true  optimum  was  less  than  that  of  a  simple  random  search. 

In  practice,  function  evaluations  (especially  those  based  on  PE)  are  time-consuming,  and  it  is  essential  to 
develop  a  strategy  which  will  have  a  rapid  rate  of  convergence,  as  well  as  a  high  likelihood  of  identifying  the 
global  optimum.  In  view  of  the  above  results,  a  two-stage  procedure  was  examined,  in  which  the  first  stage 
involved  a  random  search,  and  the  second  stage  involved  a  number  of  local  optimizations,  each  starting  from 
one  of  the  best  matches  found  during  the  search  stage.  This  two-stage  approach  had  the  advantages  of:  (1) 
the  most  thorough  sampling  of  the  parameter  space,  with  the  resulting  highest  probability  of  identifying  a 
point  in  the  region  of  the  true  optimum;  (2)  rapid  convergence  of  the  local  optimization  technique;  and  (3) 
redundancy,  when  several  different  starting  points  are  observed  to  converge  to  the  same  optimum. 

Two  local  optimization  algorithms  were  examined  for  their  effectiveness  in  this  application:  the  Powell 
direction  set  method  and  the  quasi-Newton  technique  [4].  It  was  found  that  for  independent  parameters 
(i.e.,  pij  =  0  for  ^  j),  the  direction  set  method  converged  slightly  faster,  but  when  the  parameters  were 

strongly  interdependent,  the  quasi- Newton  procedure  gave  much  faster  convergence. 

The  overall  result  of  these  experiments,  using  realizations  of  the  model  function  above,  was  the  adoption 
of  a  two-stage  strategy  involving  a  random  search  followed  by  applying  the  quasi-Newton  technique  to  opti¬ 
mize  the  best  matches  found  during  the  search. 


IV  Effectiveness  of  Two-Stage  Approach  in  Geoacoustic 
Inversion 

An  example  of  the  performance  of  this  two-stage  approach  is  shown  in  Figures  2  and  3.  Here  Orca  was  used 
to  compute  normal  mode  models  for  an  Arctic  650-m  deep  channel  at  25  Hz;  the  true  values  and  the  bounds 
on  the  parameters  to  be  optimized  are  shown  in  Table  1. 

The  two-stage  approach  was  applied  to  estimate  these  geoacoustic  parameters.  In  the  first  stage,  500 
function  evaluations  were  performed,  with  parameter  values  independently  randomly  distributed  between 
the  bounds  specified  in  Table  1.  In  the  second  stage,  the  sets  of  parameters  corresponding  to  the  10  best 
matches  (lowest  function  values)  were  used  as  starting  points  for  optimization. 
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Table  1:  Geoacoustic  parameters  (true  values  and  bounds  for  inversion)  for  the  Arctic  environment  used  for 
the  parameter  inversion  illustrated  in  Figures  2  and  3. 


Parameter  varied 

True 

Bounds 

Source-array  range  (km) 

5 

4.9 -5.1 

Water  depth  (m) 

650 

580  -  720 

Sediment  compressional  speed  (m/s) 

1800 

1620  -  1980 

Sediment  shear  speed  (m/s) 

300 

270  -  330 

Sediment  compressional  atten  (dB/wl) 

0.1 

0.07-0.13 

Sediment  shear  atten  (dB/wl) 

0.1 

0.07-0.13 

Sediment  thickness  (m) 

10 

5-15 

Basement  compressional  speed  (m/s) 

2250 

2025  -  2475 

Basement  shear  speed  (m/s) 

800 

720  -  880 

Basement  compressional  atten  (dB/wl) 

0.05 

0.03-0.07 

Basement  shear  atten  (dB/wl) 

0.20 

0.15-0.25 

Array  tilt  (‘^) 

0 

-5-5 

The  results  in  Figure  2  show  that  for  6  of  the  10  best  starting  estimates,  the  function  values  rapidly 
converged  to  zero  during  the  optimization.  For  these  cases,  essentially  perfect  matches  were  obtained  at  a 
computational  cost  of  20-30  function  evaluations  per  parameter  during  each  optimization  stage. 

The  results  in  Figure  3  show  the  evolution  of  the  12  individual  parameters  listed  in  Table  1  during 
optimization  (within  each  plot,  the  line  types  correspond  to  those  in  Figure  2).  It  is  evident  that  there 
was  a  wide  range  of  sensitivities  among  the  parameters,  with  range,  depth,  sediment  thickness,  basement 
compressional  speed,  and  array  tilt  being  most  sensitive  under  these  conditions,  and  the  others  being  almost 
invariant  as  optimization  proceeded.  For  the  sensitive  parameters,  multiple  instances  of  convergence  of  the 
different  starting  estimates  to  (approximately)  the  same  final  values  were  apparent;  these  corresponded  to 
the  convergences  of  the  function  values  to  zero  as  observed  in  Figure  2. 

Within  the  sensitive  parameters,  different  convergence  properties  were  also  observed.  For  example,  array 
tilt  converged  rapidly  to  the  true  value  in  all  cases,  while  range  converged  more  slowly  and  only  approached 
the  true  value  in  7  of  the  10  cases.  In  practice,  the  combination  of  multiple  convergences  of  the  sensitive 
parameters  to  similar  values,  combined  with  low  function  values  under  these  conditions,  should  provide  an 
effective  means  of  identifying  the  most  appropriate  estimates  for  these  parameters. 

These  results  confirm  the  effectiveness  of  the  two-stage  strategy  (i.e.,  random  search  followed  by  opti¬ 
mization  of  best  matches)  for  geoacoustic  parameter  inversion. 
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Abstract 

This  study  aims  at  the  identification  of  soTind  speed  fields  in  Mediterranean  sea  with  multisensor  analysis  by 
using  the  modal  theory.  In  Ocean  Acoustic  Tomography,  the  use  of  only  one  receiver  does  not  always  allow  an 
inversion  with  a  good  accuracy.  Therefore  a  large  vertical  array  receiver  is  used.  We  propose  here  to  realise  a 
Matched  Field  Inversion  by  using  the  beamforming  patterns.  The  principle  of  this  inversion  is  to  re-establish  the 
link  between  the  horizontal  group  velocity  and  phase  velocity  from  the  beamforming  patterns.  The  mean  sound 
speed  profile  is  reconstructed  by  using  an  adaptative  inversion  algorithm.  At  first,  we  present  the  principle  of  this 
Matched  Field  Inversion,  ctnd  then,  we  use  it  to  process  experimental  data  obtained  during  the  1994  THETIS  2 
campaign. 


1  Introduction 

Matched-field  beamforming  is  generally  used  to  determine  range/depth  localization  of  underwater  acoustic  source. 
In  Ocean  Acoustic  Tomography,  such  a  method  can  be  used  to  reconstruct  geoacoustic  properties  of  the  medium 
and  especially  the  sound  speed  profile.  In  this  study,  we  suppose  that  the  source  position  is  well  known.  We  will 
only  consider  a  range-independent  environment  and  we  will  try  to  reconstruct  the  mean  sound  speed  profile  of  the 
medium.  Generally,  the  modal  inversion  schemes  are  based  on  modal- travel  time  approaches  [1].  We  consider  here 
a  new  approach  to  process  the  inverse  problem  based  on  the  beamforming  patterns.  We  will  see  that  the  modal 
horizontal  group  velocity  and  phase  velocity  are  directly  linked  to  the  beamforming  pattern.  Using  the  group  velocity  - 
phase  velocity  velaXion^  the  inversion  algorithm  is  introduced.  A  Matched-Field  method  is  then  used  to  reconstruct  the 
mean  sound  field  by  minimizing  the  mean  square  error  between  the  experimental  and  theoritical  signal.  The  proposed 
method  allows  to  process  directly  experimental  data. 


2  Modal  approach 


Consider  a  source  at  depth  ,  which  emits  a  wideband  signal.  The  acoustic  modal  pressure  at  detph  z  from  a  distant 
source  at  range  r  can  be  expressed  as  a  sum  of  four  planes  waves  [5]: 

where  each  mode  m  is  characterized  by  the  emission  angle  6^  which  satisfies  : 


r  —  kQCOsQfn  ‘fU  Z,  _  ^ 

1  km,z  =  koSinSm 


(2) 


2.1  Group  and  phase  velocities 

The  phase  velocity  corresponding  to  the  m-th  mode  is  defined  as  vsm  =  a-nd  in  the  case  of  a  wideband  emission 


signal,  centered  around  the  angular  frequency  the  horizontal  and  vertical  group  velocities  are  respectively  defined 

.  The  group  velocity  corresponds  to  the  velocity  of  the  wavetrain 


as  v  ‘ 


h  — 
gr,m 


duj 


dk, 


m,r 


and  ± 


du! 


dk„ 


W=U)o 
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in  the  medium  ;  the  velocity  along  the  depth  for  the  vertical  group  velocity,  and  along  the  range  for  the  horizontal 
group  velocity.  Each  mode  is  propagated  through  a  waveguide  whose  thickness  hm  satisfies  c{z)  <  Vp^rn  •  The  diiferent 
types  of  propagation  can  be  classified  according  to  their  interaction  with  the  top  and  bottom  boundaries  [2]  (figure 
1).  Consider  the  ray /mode  duality  [3]  and  by  analogy  with  the  rays  theory,  the  modal  caustics  are  located  at  depth 


Figure  1:  Mode  types  according  to  their  interaction 


of  the  rays  turning  points,  A  ray,  emitted  with  an  incidence  angle  6^  (at  source  depth),  will  have  a  turning  point  at 
depth  Zm  such  that  c{zm)  =  ^^^^5  ,  where  cs  is  the  sound  velocity  at  source  depth.  Consider  the  phase  velocity  of 
the  m-mode,  we  have  : 


Cs 


=  c{Zm) 


km,r  COsOf^ 

On  a  given  receiver,  the  received  modal  angles  0^  verify  equation  (3) .  It  follows  : 


(3) 


^tp,m 


^S_Cr_ 

cose^-  cose^ 


(4) 


where  cr  is  the  sound  velocity  at  receiver  depth. 


2*2  Propagation  time 


The  modal  travel  time  is  obtained  via  the  horizontal  and  vertical  group  velocities.  At  depth  z  and  range  r,  it  can  be 
expressed  as  : 


r  ^  ^  +  e2Zs 

^gr^m  ^gr,m  (^1 ) 


with 


^^rriyZ 


\U}=UJo 


(5) 


where  Si  =  ±1  and  62  =  ±1. 

The  propagation  time  of  each  mode  depends  on  the  values  of  sj  and  62*  According  to  equation  (5),  four  propagation 
times  correspond  to  each  mode. 


2.3  Maxima  of  energy  and  group  velocity  -  phase  velocity  relation 

The  modal  positions  of  the  maxima  of  energy  of  the  received  signals  are  obtained  when  the  phase  difference  between 
two  consecutive  modes  is  a  multiple  of  27r.  This  condition  may  be  expressed  as  : 

ko[Rcos0^  +  SisinO^  {z  -j-  £2^5)]  -  ko[RcosO^_^^  SisinO^^^  {z  -f  £2^5)]  =  2p7r  (6) 


where  p  =  0, 1, 2, ...  . 

An  illustration  of  different  p  values  is  given  in  figure  2. 

Each  p  value  gives  four  modal  solutions  for  m  noted  mp(£i,  £2).  We  define  the  mean  value  Mp  of  these  four  modes  as  : 

koR{cose^^  -  cose^^^i)  =  2pTT  (7) 

By  using  equations  (4)  and  (5) ,  the  propagation  times  and  arrival  angles  of  the  received  maxima  of  energy  are  then 
written  : 


cosO^  (  \ 


c(^mp(ei,C2)) 


^mp(£i, e2)(^l5  ^2)  ^ 


(8) 


gr,mp  (e  1  ,£2)  1  2) 


(^1) 
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For  a  given  p  value,  the  four  modes  (satisfying  equation  (6)),  associated  with  the  four  maxima  of  energy,  are  neigh¬ 
bouring  and  have  approximately  the  same  physical  properties.  We  can  consider,  in  first  approximation,  that  the 
vertical  group  velocity  is  constant  in  the  vicinity  of  these  four  modes,  so  that  the  mean  value  of  the  second  term  of  the 
propagation  time  is  equal  to  zero.  The  mean  propagation  time  of  these  four  maxima  of  energy  can  be  then  written  : 


iMp  = 


^gr,Mp 


(9) 


This  travel  time  only  depends  on  the  horizontal  group  velocity.  According  now  to  equation  (4),  the  mean  received 
angles  of  the  four  maxima  of  energy  can  be  expressed  as  a  function  of  the  phase  velocity,  we  have  : 


cos9^^  — 


(10) 


From  equations  (9)  and  (10),  we  can  directly  compare  the  spatio-temporal  representation  to  the  horizontal  group 
velocity  -  phase  velocity  relation.  The  set  {£)  of  points  (^^5  cos^ . )  ^  member  of  the  horizontal  group  velocity 

-  phase  velocity  curve,  so  that,  the  different  positions  (in  the  spatio-temporal  representation)  of  the  maxima  of  energy 
of  the  received  signal  allow  to  reconstruct  the  link  between  the  horizontal  group  velocity  and  the  phase  velocity.  This 
can  be  obtained  by  using  the  following  equations  : 


^pr,n  —f 

(11) 

1 


where  n  =  1, 2, ...,  N  {N  is  the  number  of  maxima),  {in,  9n)  spatio-temporal  positions  of  the  maxima  of  energy 

at  the  receiver. 

In  other  respects,  W.  Munk  and  C.  Wunch  [3]  have  shown,  by  using  the  WKBJ  approximation,  that  the  relation 
linking  the  horizontal  group  velocity  to  the  phase  velocity  is  independent  of  the  angular  frequency  oj  and  only  depends 
on  the  sound  speed  profile.  This  group  and  phase  velocity  reconstruction  is  allowed  for  any  type  of  signal  then.  The 
result  is  totally  independent  of  the  emitted  signal  frequency. 


2  A  Simulations 

In  practice,  the  spatio-temporal  representation  is  obtained  via  a  large  vertical  array  by  beamforming.  In  this  case,  the 
cr  sound  speed  at  receiver  depth  will  be  the  mean  sound  speed  measured  on  the  whole  array. 

In  order  to  illustrate  the  group  and  phase  velocity  reconstruction  from  a  beamforming  pattern,  figure  2  represents  a 
beamforming  in  the  case  of  a  homogeneous  medium  (iso velocity  problem  c  =  1500  m/s).  We  superimpose  on  this 
picture  the  set  (X>)  of  points  i  ^  arccos  ( )J  represented  by  crosses.  The  considered  medium  is  1000  m 

^gr,m 

depth,  with  a  perfectly  reflecting  bottom.  The  source  is  placed  at  200  meters  depth,  the  receiving  array  between  300 
meters  and  500  meters  made  of  21  sensors.  The  pulsed  source  has  a  300  Hz  central  frequency,  200  Hz  bandwidth  and 
gaussian  shape.  The  propagation  length  is  50  km. 

From  the  beamforming  pattern,  if  the  mean  sound  speed  cr  along  the  vertical  array  and  the  propagation  length  are 
known,  the  spatio-temporal  positions  of  the  maxima  of  energy  allow  to  reconstruct  the  group  velocity  -  phase  velocity 
relation  (set  (i^)).  In  the  case  of  the  isovelocity  problem,  the  vertical  group  velocity  is  infinite  so  that  the  set  (^)  is  a 
member  of  the  theoretical  group  velocity  -  phase  velocity  relation  (see  figure  3,  [5]). 

In  the  general  case,  the  vertical  group  velocity  introduces  a  spatio-temporal  shift  of  the  maxima  of  energy.  The 
reconstruction  of  group  and  phase  velocities  is  obtained  by  using  equation  (11).  It  is  then  possible  to  restore  the 
horizontal  group  velocity  -  phase  velocity  relation  by  computing  the  mean  position  of  the  maxima  of  energy  for  each  p 
value.  As  an  example,  we  refer  to  the  linear  sound  speed  profile.  Figure  4  shows  the  horizontal  group  velocity  -  phase 
velocity  reconstruction  from  a  beamforming  pattern.  Figure  4  shows  clearly  that  the  set  of  mean  reconstructed  points 
{£)  is  member  to  the  theoretical  horizontal  group  velocity  -  phase  velocity  curve. 


3  Inversion  scheme 

Consider  the  horizontal  group  velocity  -  phase  velocity  relation  v<f,  =  fivgr)^  The  characteristics  of  this  curve  only 
depend  on  the  sound  speed  profile  (supposed  to  be  range-independent).  For  a  given  sound  speed  profile,  one  function  / 
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Figure  2:  Beamforming  pattern  :  case  of  the  homogeneous  medium. 


Figure  3:  Group  and  phase  velocity  relation  :  case  of  the  homogeneous  medium 
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Figure  4:  Group  and  phase  velocity  relation  :  linear  profile  case.  The  stars  represent  the  set  (5). 


exists  which  satisfies  =  fiVgr)-  Unfortunately,  the  converse  is  false.  Actually,  the  horizontal  group  velocity  -  phase 


Figure  5:  =  fi^gr)  curve  and  relation  between  the  sound  speed  and  the  thickness  of  the  waveguides 

velocity  relation  depends  on  the  different  thicknesses  of  the  modal  waveguides,  so  that,  many  sound  speed  profiles  can 
have  the  same  horizontal  group  velocity  -  phase  velocity  characteristics,  and  then  the  same  function  /  :  any  similarity 
along  the  depth  of  the  sound  speed  profile  does  not  change  the  horizontal  group  velocity  -  phase  velocity  curve.  If  the 
sound  speed  profile  is  a  monotonous  function  (it  is  the  case  in  winter  in  the  Mediterranean  sea) ,  the  thickness  of  the 
modal  waveguides  corresponds  to  the  sound  speed  profile.  In  other  cases  (presence  of  a  SOFAR),  the  thickness  of  the 
modal  waveguides  is  not  sufficient  to  reconstruct  the  sound  speed  profile. 

3*1  Algorithm 

We  developed  a  recursive  algorithm  allowing  to  reconstruct  the  relation  between  the  sound  speed  and  the  thickness  of 
the  waveguides.  The  first  step  consists  in  discretising  the  phase  velocity  z  =  0, 1,2, where  =  Co  corresponds 
to  the  minimum  value.  For  each  value  of  z,  the  algorithm  calculates  the  thickness  of  the  waveguide  (which  satisfies 
c{z)  <  up  in  order  to  obtain  the  same  horizontal  group  velocity  -  phase  velocity  relation.  The  phase  and  group 
velocities  are  calculated  with  the  KRAKEN  program. 
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3*2  Inversion  conditions 


The  accuracy  of  the  inversion  depends  on  the  knowledge  of  three  parameters  :  the  sound  speed  cr  at  receiver  depth, 
the  minimum  sound  speed  Co  and  a  particular  point  P  of  the  sound  speed  profile  (in  order  to  know  the  homothety 
ratio).  These  three  parameters  are  generally  unknown,  especially  the  cr  and  Co  parameters.  The  particular  point  P 
can  be  chosen  at  great  depth,  where  the  temperature  and  salinity  of  the  medium  are  almost  the  same  throughout  the 
year  at  a  given  depth.  It  is  then  possible  to  reduce  the  number  of  unknowns  to  two  parameters  :  cr  and  Cq.  In  order 
to  estimate  the  different  unknowns,  a  Matched  Field  Processing  method  is  used.  The  principle  consists  in  minimizing 
the  Mean  Square  Error  (MSE)  [6]  between  the  experimental  and  theoretical  time  signal,  we  have  : 

1  ^  1  f 

^(<^0  5  C/j,  P)  =  —  /  \\^exp,n{^)  ~  ^thyni^T  P)\\^dt  (12) 

where  N  is  the  number  of  sensors  on  the  array. 


3.3  Simulations 

Two  simulations  are  presented  corresponding  to  two  different  background  sound  speed  profiles.  The  results  are  pre¬ 
sented  figure  on  7  and  8.  We  consider  a  wideband  signal  centered  arround  the  central  frequency  400  Hz,  100  Hz 
bandwidth  and  gaussian  shape.  The  propagation  length  is  150  km.  The  source  is  placed  at  100  m  depth,  the  receiving 
array  between  0  and  200  m.  The  results  are  presented  in  the  case  where  the  parameters  cr  and  P  are  known. 
From  the  background  profiles  ((c),  dashed  line),  we  calculate  the  beamforming  pattern  (a)  and  the  theoritical  curve 
horizontal  group  velocity  -  phase  velocity  ((b),  solid  line).  We  now  try  to  solve  the  inverse  problem.  From  the  beam¬ 
forming  pattern,  we  directly  reconstruct  the  set  of  points  (i^))  by  using  equation  (11).  The  reconstructed  points  are 
represented  by  stars  (b).  The  inversion  algorithm  allows  to  reconstruct  the  sound  speed  profile  ((c),  solid  line),  A 
comparison  of  both  simulations  shows  that  the  inversion  accuracy  is  very  satisfying.  In  the  case  where  parameters  c^, 
Cr  and  P  are  known,  the  calculation  cost  is  low  (less  than  5  minutes).  In  other  cases,  we  use  the  method  described  in 
the  previous  section.  The  calculation  cost  is  far  bigger  then,  but  the  inversion  accuracy  is  the  same  (in  simulations). 
This  inversion  method  allows  to  reconstruct  the  sound  speed  profile  without  any  a  priori  information  and  can  be 
directly  applied  to  experimental  data. 

3.4  Experimental  data  processing 

We  applied  this  inversion  method  to  experimental  data  obtained  during  the  1994  THETIS  2  campaign.  We  consider 
here  the  propagation  between  a  mooring  (W5)  situated  close  to  the  Corse  coast  and  a  reception  array  made  of  16 
sensors,  situated  close  to  the  French  coast  (Cap  Ferrat).  We  processed  here  only  one  record  dated  the  28th  of  September 
1994  at  llh28  UTC.  The  array  characteristics  are  as  follows  :  the  first  hydrophon  is  situated  at  depth  145,8  m  and 
the  last  one  at  184.2  m,  the  14  others  are  vertically  distributed  between  these  two  depths.  At  the  record  date,  the 
source  depth  (W5)  is  125.5  m,  and  the  propagation  length,  calculated  with  the  Andoyer’s  formula,  is  153.427  km. 
The  emited  signal  is  a  series  of  15  or  30  sequences  of  5.11  seconds  BPSK.  In  this  propagation  configuration,  we  try 
to  determine  the  mean  sound  speed  field  between  W5  and  the  vertical  array.  Unfortunately,  the  used  vertical  array 
is  too  short  :  we  have  no  information  between  0  m  and  145  m,  and  especially  at  the  SOFAR  depth  (around  100  m 
depth),  so  that  it  is  almost  impossible  to  reconstruct  the  sound  speed  between  the  surface  and  the  SOFAR  depth. 

By  using  the  inversion  scheme  and  by  minimizing  the  Mean  Square  Error  defined  by  equation  (12),  it  has  been  possible 
to  reconstruct  the  sound  speed  field  for  only  deep  water  (from  100  m  to  2000  m).  The  first  100  m  were  reconstructed 
by  using  experimental  measurements  obtained  in  October  1994,  50  km  from  the  array.  The  beamforming  pattern  and 
the  inversion  result  are  presented  in  figure  9.  The  recorded  and  predicted  arrival  patterns  are  compared  (for  the  central 
receiver). 
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4  Source  or  receiver  detph  estimation 


It  is  possible  to  estimate  the  depth  of  one  sensor  (transmitter  or  receiver)  by  using  a  group  of  four  maxima  of  energy. 
The  group  velocity  corresponding  to  a  maximum  of  energy  is  given  by  : 


£l(z  +  £2Z,) 


(13) 


Using  the  approximation  that  the  vertical  group  velocity  is  approximatively  the  same  around  the  four  maxima,  we 
have  : 

t;^r(£l,£2)  ~  Vpr  (14) 

It  is  then  possible  to  express  the  depth  of  an  instrument  with  respect  to  another.  By  using  the  four  propagation  times 
fhe  source  depth  according  to  the  receiver  depth  is  : 

and  we  can  express  the  receiver  depth  according  to  the  source  depth  by  : 

,  ^  [^(1, 1)  -  t(-i,  1)]  +  [^(1,  -1)  -  f(-i,  -1)] 


In  order  to  apply  this  processing  to  experimental  data,  we  do  not  use  the  absolute  propagation  times,  but  only  the 
relative  propagation  times.  The  result  is  totally  independent  of  the  propagation  length  r.  If  the  propagation  times 
are  known  to  within  1  ms,  the  estimation  error  is  less  than  1  %.  Such  a  method  may  be  used  to  estimate  the  source 
depth  or  the  receiver  depth  (for  example  when  a  synthetic  array  is  used)  [5], 


Figure  6:  representation  of  a  group  of  four  maxima  of  energy 


5  Conclusion 

In  this  study,  we  have  seen  that  the  beamforming  pattern  allows  to  reconstruct  the  horizontal  group  velocity  -phase 
velocity  relaXion.  The  main  application  of  such  a  processing  is  the  estimation  of  the  sound  speed  profile.  The  Matched 
Field  Processing  method,  presented  in  chapter  3.2,  is  then  used.  This  inversion  method  may  be  directly  used  to  process 
experimental  data  and  allows  to  reconstruct  the  sound  speed  profile  without  any  a  priori  information.  The  inversion 
accuracy  depends  on  the  number  of  sensors  and  the  length  of  the  vertical  array.  In  the  Mediterranean  basin,  we  can 
hope  a  very  good  inversion  accuracy  by  using  a  200  m  high  array  from  the  surface  to  200  m  depth.  Unfortunately,  it 
was  not  the  case  in  our  experimetal  data.  Another  application  is  the  source  or  receiver  depth  estimation  presented  in 
chapter  4.  In  the  future,  we  will  try  to  estimate  geoproperties  of  the  bottom  by  using  the  horizontal  group  velocity  - 
phase  velocity  reconstruction  and  the  inversion  by  Matched  Field  Processing. 
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Figure  7:  (top)  beamforming  pattern  obtained  with  the  background  profile.  (Down  left)  Group  and  phase  velocities 
reconstruction  (the  dashed  line  corresponds  to  the  line  ^  Vgr)^  (Down  right)  Inversion  result. 
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Figure  8;  (top)  beamforming  pattern  obtained  with  the  background  profile.  (Down  left)  Group  and  phase  velocities 
reconstruction  (the  dashed  line  corresponds  to  the  line  =  Vgr).  (Down  right)  Inversion  result. 
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Figure  9:  (top)  beamforming  pattern  -  (down  left)  reconstructed  sound  speed  profile  -  (down  right)  experimental  signal 
(solid  line),  reconstructed  signal  (dash  dot) 
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Abstract 

For  over  a  decade,  advanced  signal  processing  techniques  based  on  Matched-Field,  Matched-Mode 
or  Model-Based  Processing  (MFP/MMP/MBP)  concepts  have  been  developed  for  improving  the  ca¬ 
pability  of  passive  and  active  sonar  systems.  Studies  carried  out  at  various  laboratories  have  shown 
that  MFP/MMP/MBP  techniques  can  provide  increased  performance  over  conventional  beamforming 
schemes  for  detecting  and  localizing  underwater  targets.  At  the  Esquimalt  Defence  Research  Detach¬ 
ment,  a  second  generation  of  MFP/MMP  based  software  called  the  Advanced  Acoustic  Signal  Processor 
(AASP)  is  currently  being  developed  in  collaboration  with  contractors  and  researchers  at  universities. 
The  basic  components  of  the  AASP  system  consist  of  data  preconditioning,  an  environmental  database, 
a  propagation  model,  a  matcher  and  a  tracker.  An  inversion  algorithm  is  also  included  for  updating  the 
environmental  parameters.  Hydrophone  signals  received  from  a  surveillance  or  tactical  array  are  corre¬ 
lated  with  replica  fields  generated  using  a  propagation  model  to  cover  the  search  region.  This  matching 
is  carried  out  for  many  candidate  target  locations  to  form  an  ambiguity  surface  whose  values  reflect 
the  likelihood  that  the  source  is  at  the  positions  in  the  search  region.  A  sequence  of  such  ambiguity 
surfaces  in  time  is  processed  by  a  detect-after-track  algorithm.  Uncertainties  in  the  knowledge  of  the 
environmental  parameters  can  degrade  the  matching  between  the  measurements  and  the  replicas:  these 
can  be  mitigated  by  the  tracker  or  reduced  by  using  geoacoustic  inversion  methods  to  focus  the  ambiguity 
surfaces.  The  AASP  is  being  designed  to  run  on  shared-memory  and  distributed  processors  connected  by 
a  local  area  network  and  will  incorporate  the  lessons  learned  from  the  previous  generation  of  MFP/MMP 
systems  into  an  end-to-end  research  tool  for  the  development,  testing  and  evaluation  of  individual  pro¬ 
cessor  components.  In  this  paper,  the  nature  of  the  processor  components  and  their  integration  into  the 
AASP  will  be  discussed. 


Introduction 

The  objective  in  developing  the  Advanced  Acoustic  Signal  Processor  (AASP)  is  to  produce  portable  software 
for  research  and  development  of  Matched-Field  Processing  (MFP).  This  is  to  be  achieved  through  a  collab¬ 
orative  effort  that  involves  the  Department  of  National  Defence  (DND),  the  private  sector  and  universities. 
The  funding  from  DND  is  subsidized  by  reduced  rates  from  the  prime  contractor  and  makes  use  of  expertise 
existing  among  the  subcontractors.  Economy  is  also  achieved  through  reuse  of  existing  code  held  at  the 
Esquimalt  Defence  Research  Detachment,  a  division  of  Defence  Research  Establishment  Atlantic  (DREA) 
and  at  DREA. 
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Previously  our  Matched-Field  Processing  (MFP)  has  been  carried  out  in  stages  with  significant  operator 
intervention.  The  development  of  an  end-to-end  AASP  software  package  will  enable  analysis  of  data  with  less 
skilled  manpower.  This  is  possible  on  account  of  the  more  advanced  state  of  MFP  and  software  development 
environments  and  is  a  natural  step  in  the  development  of  a  prototype  on  the  route  to  an  operational  system. 

The  software  is  being  developed  in  a  modular  form  to  enable  subsequent  enhancements  to  a  variety  of 
applications.  The  first  priority,  that  is  anticipated  will  take  two  years  to  complete,  is  the  development  of 
surveillance  oriented  software  for  a  fixed  array.  It  is  planned  that  subsequent  enhancements  to  the  software 
for  data  from  sonobuoys,  geobuoys  and  towed  arrays  will  follow. 

The  research  on  which  the  AASP  is  based  is  documented  in  numerous  papers  and  is  broadly  supported 
by  the  MFP  literature.  Recently,  the  emphasis  at  EDRD  has  been  to  investigate  those  aspects  necessary  for 
the  development  of  an  operational  system.  Results  from  data  analysis  obtained  during 

the  past  two  years  indicate  that  it  is  now  appropriate  to  develop  a  prototype.  Such  a  prototype  will  allow 
the  timely  evaluation  and  demonstration  of  end-to-end  MFP  software  for  R&D  that  is  the  last  step  before 
developing  an  experimental  development  model. 

The  remainder  of  this  paper  describes  the  environmental  database,  propagation  modelling,  matching, 
tracking  and  inversion  components,  and  gives  the  reasoning  behind  the  choices  made.  This  is  followed  by  a 
description  of  the  implementation  of  the  rapid  prototype  and  its  capabilities.  Some  details  of  the  propagation 
modelling  will  be  included  while  such  details  for  other  components  may  be  found  in  the  references  cited. 


System  components 

Overview 

MFP  consists  of  matching  acoustic  data  with  predictions  of  the  data  for  all  possible  source  positions  to  find 
the  matching  source  range,  depth  and  bearing.  This  leads  to  improved  array  gain  and  a  better  determination 
of  the  source  position  than  plane  wave,  bearings-only,  beamforming.  In  the  first  step  in  MFP,  acoustic 
data  is  windowed  and  Fourier  transformed  prior  to  MFP  in  a  similar  manner  to  that  used  in  plane  wave 
beamforming.  Replica  vectors,  predictions  of  the  received  data,  for  all  source  positions  in  the  search  region  are 
then  calculated  from  the  propagation  model  using  an  estimate  of  the  environment  drawn  from  the  database. 
Each  of  these  replicas  is  compared  to  a  measured  data  vector  to  form  points  on  an  ambiguity  ’surface’.  A 
sequence  of  data  vectors  results  in  a  sequence  of  ambiguity  surfaces.  In  our  implementation,  linear  or  circular 
tracks  are  formed  through  the  largest  peaks  on  the  sequence  of  ambiguity  surfaces  to  determine  the  track 
with  the  highest  signal-to-noise  ratio  and  therefore  the  most  likely  track. 

MFP-based  processing  on  sparse  environmental  data  has  demonstrated  detection  of  weak  sources  and  the 
determination  of  their  range,  depth  and  bearing.  Further  improvement  in  performance  can  be  expected  when 
a  better  model  of  the  environment  is  available.  MFP  inversion  for  estimating  the  environmental  properties 
can  be  carried  out  on  a  continuing  basis  for  this  purpose. 

Environmental  database 

We  have  employed  files  to  store  environmental  data  for  propagation  modelling  in  the  past  with  the  user 
entering  data  for  each  application.  This  has  been  a  laborious  process  that  has  been  repeated  for  each 
application  with  little  use  being  made  of  existing  databases  because  of  the  large  initial  cost.  A  more  efficient 
method  of  storing  and  retrieving  environmental  data  is  desirable  now  that  repeated  database  use  is  likely.  A 
high  performance  object-relational  spatial  database  management  system  called  ILLUSTRA  has  been  chosen 
for  the  development  of  such  a  database.  Because  ILLUSTRA  is  a  spatial  database,  it  is  efficient  in  two  and 
three  dimensional  searches  of  environmental  data  required  by  the  AASP.  Subsequent  development  to  provide 
a  database  for  all  the  MFP  data  files  may  also  be  carried  out  with  ILLUSTRA. 

Initially,  the  Environmental  DataBase  (EDB)  is  to  include  both  environmental  data  and  geometric  data 
pertaining  to  the  array  position.  Included  in  the  environmental  data  are  bathymetry  as  well  as  sound  speed(s) 
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and  attenuations  in  the  water,  bottom  and  ice.  Additional  parameters  such  as  boundary  roughness  could  be 
added  later. 

The  EDB  has  a  layered  modular  design  [1,  2].  There  are  three  modules,  a  data  visualization  module, 
an  interpolation  module  and  the  ILLUSTRA  module  itself.  The  interpolation  module  carries  out  the  two 
and  three  dimensional  interpolation  to  provide  the  data  files  required  by  the  propagation  model.  Data 
visualization  is  used  to  check  on  features  of  the  data  stored  in  the  database  or  for  presentation  of  the  data. 
This  is  particularly  useful  when  new  data  has  been  entered  into  the  database  and  some  check  of  the  success 
of  the  process  is  required. 


Propagation  model 

Numerical  PE  predictions  are  routinely  used  to  model  underwater  sound  propagation.  The  main  rationale 
for  this  is  that  accurate  full-wave  solutions  to  the  PE  can  be  computed  efficiently  using  marching  algorithms 
for  both  depth-  and  range-dependent  inhomogeneous  media  [3].  The  development  of  the  PE  method  has 
reached  the  point  where  finite-difference  implementations  derived  from  Pade  series  expansions  can  provide 
accurate  solutions  to  one-way  wave  propagation  for  realistic  geoacoustic  conditions,  including  variable-depth 
bathymetry.  Moreover,  with  the  introduction  of  both  exact  and  approximate  PE  procedures  for  handling 
elastic  media,  the  physics  of  shear  wave  propagation  can  be  incorporated  in  a  straightforward  way. 

In  two  dimensions,  range  r  and  depth  z  with  positive  down,  let  the  region  0  <  z  <  zi,  below  the 
ocean  surface  {z  =  0)  be  occupied  by  a  fluid  with  density  p(r,0),  compressional  sound  speed  c{r,z)  and 
compressional  absorption  a(r,  z).  For  z  >  z^,  these  material  properties  are  assumed  to  have  the  constant 
values  pby  Cb  and  respectively.  In  the  case  of  a  solid  medium,  the  shear  speed  Cg  and  shear  absorption 
as  must  also  be  specified.  It  is  convenient  to  define  N  =  n(l  -h  ia)  where  n  =  co/c  is  the  usual  acoustic 
refractive  index,  and  A^o  =  t<^/co  is  a  reference  wavenumber. 

For  r  >  0,  the  outgoing  component  of  the  spatial  pressure  field  p(r,z)  =  '0(r,z)  exp(iA^r)/*;/r  produced 
by  a  harmonic  exp{—iujt)  point  source  located  at  (r,  z)  =  (0,Zs)  can  be  recovered  from  the  step-by-step 
numerical  solution  of  the  higher-order  Pade  PE 
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where  the  Pade  coefficients  aj^j,  bj^j  and  the  operator  are  defined  by  [4] 
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Equation  (1)  can  be  solved  for  the  field  at  r  +  Ar  in  terms  of  the  field  at  r  as  a  sequence  of  J  systems  of 
equations,  where  the  yth  system  is  given  by 
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and  where  we  have  set  cfj  =  bjj  ±  ^ikoArUj^j ,  If  —  l)  ^(r,  z)  is  evaluated  in  the  heterogeneous 
approximation  [5]  at  each  point  z  on  the  depth  grid,  namely 
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where  p±  and  po  denote  the  density  combinations 
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then  Eq.  (3)  becomes  a  tri-diagonal  system  of  equations  that  can  be  efficiently  solved  by  double  recursion. 
A  simple  analytic  starter  is  used  to  specify  an  initial  field  at  r  =  0.  The  radiation  condition  at  z  oo  is 
approximated  by  appending  an  absorbing  layer  to  the  computational  grid  and  applying  a  pressure-release 
boundary  condition  at  the  base  of  this  layer. 

Two  aspects  of  the  above  implementation  deserve  comment.  First,  it  is  convenient  in  some  applications 
to  replace  the  real  coefficients  ajj  and  bjj  in  Eq.  (2)  with  complex  ones.  For  j  <  8,  suitable  values  have  been 
determined  numerically  and  tabulated  [6].  Second,  although  a  PE  formalism  exists  for  solving  the  elastic 
PE  [3,  6,  7],  it  is  more  efficient  to  make  use  of  an  equivalent  fluid  approximation  developed  by  Zhang  and 
Tindle  [8].  In  this  case,  the  (real)  density  p  of  an  elastic  medium  is  simply  replaced  by  a  complex  density  p* 
that  is  determined  from  the  shear  properties  of  the  medium  according  to 


p*{r,z)  ^p{r,z) 
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where  Ng  =  (co/cs)(l  +  ia^). 


Matching 

The  matcher  is  a  critical  component  in  the  determination  of  the  robustness  of  the  MFP  software.  As 
with  all  components,  the  simplest  matcher  that  provides  the  detection  and  localization  capability  required 
was  selected.  For  low  Signal-to-Noise  ratios  (SNR’s)  in  spatially  uncorrelated  noise  whose  amplitudes  are 
Gaussian  distributed,  the  Bartlett  matcher  represents  an  optimum  matcher.  Since  this  is  simply  an  inner 
product  of  the  Fourier  transformed  data  vector  and  the  replica  vectors,  it  is  also  easy  to  implement.  Such 
a  choice  implies  that  the  data  needs  to  be  prewhitened  for  the  prevailing  noise.  It  would  also  be  desirable 
to  remove  impulses  in  the  noise  in  order  to  satisfy  the  assumption  that  the  noise  amphtudes  are  Gaussian 
distributed.  Lastly,  the  replicas  are  normalized  to  have  unit  norm  to  ensure  that  in  spatially  white  noise  the 
ambiguity  surfaces  are  flat. 

Since  the  computational  load  of  matching  would  limit  the  applicability  of  MFP,  schemes  have  been 
investigated  for  reducing  the  load  without  incurring  an  unacceptable  performance  degradation.  The  per¬ 
formance  enhancements  to  be  expected  at  high  signal  to  noise  ratios  are  substantial  [9].  At  low  SNR’s  the 
computational  load  can  be  reduced  by  an  order  of  magnitude  at  a  cost  of  a  1-dB  SNR  penalty  [10,  11] 

Tracking 

At  high  SNR’s,  and  when  the  environmental  parameters  are  sufficiently  well  known,  the  largest  peak  in  an 
ambiguity  surface  will  occur  at  the  true  source  position  with  sidelobes  at  a  lower  level.  The  position  of  these 
sidelobes  will  depend  on  the  source  position  and  the  scenario.  At  low  SNR’s,  or  when  the  environment  is 
poorly  known,  the  largest  peak  may  not  occur  at  the  true  source  position.  On  average  the  largest  peak 
will  be  found  most  frequently  at  the  source  position.  If  the  source  is  moving,  the  sidelobes  will  vary  in 
their  position  relative  to  the  source  position.  Thus  integrating  the  Bartlett  statistic  on  successive  ambiguity 
surfaces  along  all  possible  source  tracks  would  allow  one  to  recover  the  source  track,  as  it  would  have  the 
highest  integral.  In  practice  this  takes  too  long  because  of  the  enormous  number  of  possible  source  tracks. 
Our  solution  is  to  integrate  along  linear  tracks  which  pass  through  the  largest  peaks  on  pairs  of  ambiguity 
surfaces  [12,  13,  14].  This  reduces  the  computation  time  about  seven  orders  of  magnitude  for  search  regions 
of  interest.  Our  efficient  tracking  algorithm  also  determines  the  track  segment  over  which  a  source  appears 
to  be  present.  The  present  version  of  the  tracking  algorithms  search  for  linear  or  circular  track  segments. 
These  are  then  concatinated  to  determine  the  overall  track  shape.  This  efficient  algorithm  incurs  only  a 
small  performance  loss  compared  to  an  exhaustive  search  and  yet  can  be  implemented  in  a  real  time  system. 


Inversion 

Inversion  consists  of  finding  the  set  of  environmental  and  geometric  parameters  that  provide  a  replica  that 
best  matches  the  data.  As  with  tracking,  this  can  be  computationally  prohibitive.  Restricting  the  number 


of  environmental  parameters  to  those  that  most  affect  the  match  reduces  the  computational  load.  Synthetic 
anealing  and  genetic  algorithms  have  been  widely  reported  on  for  the  purpose.  Recently,  a  two-step  method 
has  been  investigated.  In  the  first  step,  candidates  for  the  optimum  parameter  set  are  found  through  a 
random  search.  This  is  followed  by  a  gradient  search  to  find  optimum  values.  This  two-step  algorithm 
has  been  compared  favorably  to  the  synthetic  anealing  and  genetic  algorithms  in  that  it  finds  the  optimum 
parameter  set  more  efficiently  and  more  consistently  in  simulations  [15].  The  gradient  search  also  has  the 
advantage  of  requiring  no  user  specified  parameters.  Synthetic  anealing,  genetic  and  two-stage  searches  will 
be  investigated  more  fully  with  the  prototype. 


Implementation  and  system  capabilities 

Implementation 

In  order  to  minimize  programming  cost  and  yet  provide  a  user-friendly  software  package,  the  software  was 
and  is  being  written  in  Research  Systems  Interactive  Data  Language  (IDL).  This  allows  efficient  development 
of  widget-driven  software  that  minimizes  the  need  for  an  extensive  users  guide  for  a  novice  operator.  This 
software  approach  also  minimizes  data  entry  and  makes  use  of  a  point  and  click  approach  for  efficiency  of 
use.  IDL  code  may  be  executed  on  a  wide  range  of  single  CPU  platforms  and  can  provide  a  good  graphical 
output  at  a  reasonable  programming  cost.  Furthermore  code  written  in  IDL  is  more  easily  maintained  by  a 
mixed  group  of  defence  scientists  and  contractors  at  separate  locations.  For  the  protyping  of  any  component, 
it  was  also  considered  too  time  consuming  to  write  code  in  a  language  such  as  FORTRAN.  Instead,  it  was 
decided  to  implement  all  code  directly  in  IDL.  Where  necessary  existing  components  in  other  languages 
could  be  accommodated  and  resources  permitting  other  more  computationally  efficient  codes  would  be  used 
for  the  kernels  in  an  experimental  demonstration  model.  Our  code  has  been  tested  under  Unix  on  HP,  Sun 
and  SGI  work  stations  as  well  as  under  Windows  on  PCs  and  MacOs  on  the  Power  Macintosh. 


System  capabilities 

At  present,  a  rapid  prototype  exists  that  allows  simulation  and  analysis  of  data 

for  a  vertical  array.  Simulation  includes  the  addition  of  spatially  white  noise  to  the  signal  obtained  from 
the  propagation  model  while  analysis  includes  matching,  tracking  and  inversion  for  environmental  properties. 
The  propagation  model  assumes  a  range-independent  environment  and  the  propagation  modelling  is  N  by 
2-D.  The  matching  is  exhaustive  while  the  tracking  assumes  a  source  moving  along  a  straight  line  or  circle 
in  three  dimensions.  Enhancements  planned  include  3-D  propagation  modelling,  the  replacement  of  the  file- 
based  database  with  the  ILLUSTRA-based  database,  efficient  matching,  and  the  development  of  a  version 
for  a  multi-CPU  platform.  The  intention  is  that  the  prototype  enable  the  efficient  evaluation  of  these  and 
other  enhancements. 


Concluding  remarks 

Experience  to  date  with  the  preliminary  version  of  the  software  shows  that  it  has  proven  portable,  to  have  a 
shorter  learning  curve  for  its  use,  to  have  relatively  trouble-free  execution  and  is  easy  to  modify.  It  has  the 
look  and  feel  of  recently  developed  userfriendly  code.  As  expected  it  executes  more  slowly  than  code  written 
in  FORTRAN.  This  is  an  acceptable  tradeoff  for  a  development  prototype.  The  present  version  took  about 
one  man  year  to  produce  and  contains  many  of  the  features  planned  for  the  final  system.  Three  dimensional 
propagation  modelling,  a  true  database  and  a  multiprocessor  capability  are  the  principle  developments 
planned  in  the  next  phase. 
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